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THREE-DIMENSIONAL ANALYSIS OF 
SURFACE CRACK-HERTZIAN STRESS FIELD INTERACTION 


R. Ballarini and Y. Hsu 
Dept, of Civil Engineering 
Case Western Reserve University 
Cleveland Ohio 

SUMMARY 

This thesis presents the results of a stress intensity factor analysis of 
semicircular surface cracks in the inner raceway of an engine bearing. The loading 
consists of a moving spherical Hertzian contact load and an axial stress due to 
rotation and shrink fit. 

A three dimensional linear elastic Boundary Element Method code was 
developed to perform the stress analysis. The element library includes linear and 
quadratic isoparametric surface elements. Singular quarter point elements were 
employed to capture the square root displacement variation and the inverse square 
root stress singularity along the crack front. The program also possesses the 
capability to separate the whole domain into two subregions. This procedure 
enables one to solve non-symmetric fracture mechanics problems without having to 
separate the crack surfaces a priori. 

A wide range of configuration parameters was investigated. The ratio of 
crack depth to bearing thickness was varied from one-sixtieth to one-fifth for 
several different locations of the Hertzian load. The stress intensity factors for 


v 


several crack inclinations were also investigated. 


The results demonstrate the efficiency and accuracy of the Boundary 
Element Method. Moreover, the results can provide the basis for crack growth 
calculations and fatigue life prediction. 



CHAPTER ONE 


INTRODUCTION 


Surface cracks commonly occur in machine and structural components. An 
example of such a component is a rotating engine bearing subjected to rolling 
contact. Under high speed rotation and cyclic contact loading, the surface crack 
initiating at the raceway of the bearing might propagate and lead to catastrophic 
failure. Raceway fracture is a totally unacceptable failure mechanism because it 
may cause serious damage to engine operations and consequently produce 
catastrophic engine failure. An accurate crack stress analysis of the surface-cracked 
component is essential in order to make a reliable prediction of fatigue life. 
However, due to the complexities of the nature of the surface crack problem, 
mathematical closed form solutions are not possible, and a numerical analysis or an 
experimental approach must be used to determine the stress intensity factors for 
surface cracks under different types of loading. The Boundary Element Method is 
an efficient and accurate tool for fracture mechanics analyses if singular elements 
and multi-domain crack modeling are employed. This method is used in this 
research. 

Several factors will affect the growth of surface cracks in a rotational 
bearing under rolling contact loads. These include the geomety and inclination of 
the crack, the tensile hoop stress due to rotation and shrink fit, the moving Hertzian 
load, the pressure of the lubricant seeping into the crack, the shear stress on the 
raceway surface due to the sliding contact, and friction along the crack surfaces. A 
significant amount of research has been conducted aimed at gaining a better 
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understanding of the effects of each of these factors. While the surface crack is a 
three-dimensional problem, most of the analyses which appear in the literature are 
two-dimensional. These include the work of Way [1], which considers the effect 
of the lubricant, Fleming and Suh [2,3], which considers the effects of surface 
friction, Rosenfield [4], which considers the effects of crack surface friction, and 
Clark [5], which considers the effects of tensile hoop stresses. A recent paper by 
Mendelson and Ghosn [6] presents the results of fatigue life predictions of a 
propagating surface crack subjected to tensile hoop stresses and cyclic Hertzian 
contact loadings. Using a modified Forman-type crack propagation law they 
predicted the fatigue life of a typical bearing and compared their results with 
experimentally observed fatigue lives. Their predictions were conservative by a 
factor of 12. However, they demonstrated that the crack driving force in such 
problems is the alternating mixed-mode loading that occurs with each passage of the 
roller. Based on these results, the present research was aimed at quantifying the 
three-dimensional effects of the problem. Three dimensional analyses of surface 
cracks as applied to contact fatigue were recendy performed by Murakami [7]. 
However, in his analysis the tensile hoop stresses were ignored. The model 
proposed in the present research neglects some of the factors mentioned previously. 
It is assumed that lubrication renders surface sliding friction negligible. The 
pressure on the crack surfaces which may arise from the lubricant seeping into the 
crack is ignored since the Hertzian loading moves past the crack very fast, thus the 
viscosity, compressibility, and inertia of the oil will prevent pressurization of the 
crack surfaces [8]. Moreover, since the radius of the Hertzian contact area is 
smaller than the surface length of the crack, the crack mouth will not be completely 
covered and the oil is allowed to squeeze out of the crack. The friction between the 
crack surfaces is neglected, since it tends to increase the resistance to crack growth. 
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This will lead to a conservative prediction of fatigue life. 


Thus, the only factors assumed to be important in the present model are the 
mechanical loads arising from the Hertzian contact, rotation, and shrink-fit. The 
remaining five chapters of this paper are organized as follows: 

In Chapter two the boundary integral equation is derived. The equation is 
reduced to a system of algebraic equations, and a procedure is described which 
treats the singularities which appear in the kernels. An algorithm for multi-domain 
analyses is also presented. 

Chapter three is a brief review of linear elastic fracture mechanics. The 
formation of quarter point elements and traction singular elements as well as the 
displacement correlation method for calculating the stress intensity factors are also 
discussed. Several verification problems follow in Chapter four to elucidate the 
accuracy and efficiency of the Boundary Element Method for solving three 
dimensional linear elastic solid mechanics problems including crack stress analysis. 

In Chapter five, the spherical Hertzian stress distribution and the hoop 
tensile stress due to the rotation and shrink fit of the inner raceway of the engine 
bearing are calculated. Results are presented for a wide range of configuration 
parameters. These include several different locations of the Hertzian load, different 
inclinations of the crack surface, several ratios of the crack depth to the raceway 
thickness and different intensities of the Hertzian load. A large number of the stress 
intensity factor versus these factors are presented. 
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The final chapter presents a discussion of the results and recommendations for 
future research. 
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CHAPTER TWO 


BOUNDARY ELEMENT METHOD 


This chapter reviews the development of the Boundary Element Method and 
presents a detailed derivation of the boundary integral equation. The procedures of 
numerical implementation and a multi-domain technique of the Boundary Element 
Method are also illustrated. 

2.1 Introduction 

The boundary integral equation was first derived explicitly by Rizzo [9], 
who reduced two dimensional isotropic elastostatics problems to an integral 
equation by using the Betti-Somigliana formula. The equation was then discretized 
into line segments along the boundary over which constant displacement and 
traction distributions are assumed to solve the elastic problem numerically. 
Although much of the mathematical theory of Rizzo's formula can be traced to 
Kupradze [10], and several papers about integral equation methods such as Jawson 
[11] and Sy mm [12], his work presented a clear form of the integral equation 
relating the boundary traction and displacement which is now commonly associated 
with the term " Boundary Integral Equation". Following Rizzo’s work, Cruse [13] 
extended the Boundary Element Method to three dimensional isotropic elastic 
problems by using triangular elements with linear variations. The accuracy of the 
Boundary Element Method was later improved by Lachat and Watson [14], and 
Rizzo and Shippy [15] using higher order isoparametric elements. The application 
of the Boundary Element Method to fracture mechanics was carried out by Cruse 
[13], [16] by modeling the crack as an open notch. The results obtained using this 
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approach tend to be inaccurate because as the cracks surfaces are moved close to 
each other the system of equations becomes singular. This problem was remedied 
by Blandford et al. [17], who used two subregions to model the crack. The 
quarter-point technique developed for the Finite Element Method by Barsoum [ 1 8] 
was adopted and modified by Cruse and Wilson [19] to capture the square root 
singularity predicted by linear elastic fracture mechanics. 

The Boundary Element Method developed in this paper utilizes four kinds 
of isoparametric surface elements: a three-node linear triangular element, a four- 
node linear quadrilateral element, a six-node quadratic triangular element, and eight- 
node quadratic quadrilateral element. A library of Gaussian integration quadrature 
is installed in a subroutine which can be used to accomplish the numerical 
integration. The quarter point element and the traction singular element are used in 
the program to represent the crack tip singularity. The multi-domain technique is 
also applied to model the topology. 


2.2 Derivation of Boundary Integral Equation 


The mathematical fundations of the boundary integral equation are based on 
the Kelvin solution and Betti's reciprocal theorem [20]. 

Let P and Q be two arbitrary points in an infinte elastic body as shown in 
Fig 2.1. A unit concentrated load acting at point P in i direction is defined as 
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Figure 2. 1 Domain and Boundary Geometry 
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fi (P,Q)= 5(P,Q) ei 


where is the unit vector in i direction and 5(P,Q) is the Dirac delta function which 
is defined as a function that is equal to zero for Q does not coincide with P and 
becomes infinite when P=Q in such manner that 


J 5(P,Q) dfl(Q) = 1 
n 


for any point P which lies in domain £2. The displacement in direction j at point Q 
due to a unit concentrated load applied in direction i at point P in an infinte linear 
elastic body is given by Kelvin's solution [21] 


u. 

j 


= U M (P ,Q) e. 


( 2 . 1 ) 


where 


U..(P,Q)* 


1 

16rcp.(l-v)r 


[ ( 3-4v)5..+r, i r, j ] 


( 2 . 2 ) 


and e^ is the component of the unit base vector in direction i, jj. is shear modulus, v 
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is Poisson’s ratio, and r is the distance between point P and point Q. Note that 


r= 


. 1/2 


(2.3) 


and 




r 


r = &_ = 1(5, 

’ij dx.dx. r K % 


r,.r,. ) 

’l 'j ' 


(2.4) 


where 5jj is the Kronecker delta which is equal to 1 when i = j and 0 when i * j and 
comma i denotes partial differentiation with respect to direction xj. Cartesian tensor 
notation is used hereinafter with all the subscripts ranging from 1 to 3 and the 
convention that repeated indices are summed is employed. The stress-strain relation 
for an isotropic linear elastic material is [22] 



and the strain -displacement relation is 


e jk 2 V + 



( 2 . 6 ) 


The stress-displacement relationship for an isotropic elastic material are obtained by 
substituting Eq.2.6 into Eq.2.5, i.e. 




+ M u yk + 




(2.7) 


where X = ^ ^ is ^ oun §’ s Modulus 


Substituting Eq.2.1 into Eq.2.7, the stress field can be obtained as 


D ijk( P 'Q) = U^, m (P ,Q) + n( U^(P,Q) 


+ U^CP-Q)) 


-l 

8 tc ( 1 — v )!- 2 



( 2 . 8 ) 
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where Djj^ (P,Q) is interpreted as the stress component at point Q due to a unit 
load in the i direction at point P. If point Q is put on the boundary T of a finite 
body with domain Cl cut out from the infinite body.as shown in Fig.2.1, the 
tractions at point Q on the surface can be determined as 


t. 

j 


= Vk 


(2.9) 


where n^ is the component of the outward normal to the surface at point Q. 
Substitution of Eq.2.8 into Eq.2.9 leads to 


T yCP.Q> " VRQK(Q> 


8tc (l-v)r 2 


(l-2v)(n j r,.-n i r, j ) - n k r, k f(l-2v)5. j + 3r,.rJ|(2.10) 


where Tij(P,Q) is the traction at point Q in direction j on the surface with outward 
normal n^ due to a unit load in direction i at point P. The free body cut out from the 
infinite body forms an equilibrium state subjected to the concentrated unit force 
fj(P.Q) and the boundary tractions Tjj(P,Q)e^ with the corresponding boundary 
displacement Uij(P,Q)ep Betti's reciprocal theorem can now be applied to derive 
the boundary integral equation. Suppose that there are two generalized force 
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systems. The first system includes body forces bj, surface tractions tj and 
displacements uj , and the second one consists of body forces bj*, boundary 
tractions tj and displacements Uj . If these two systems act simultaneously on a 
linear elastic body with domain Q enveloped by the boundary surface F, Betti's 
reciprocal theorem states that [23] 


Jt. u*dT + J b.u* dO = Jt* u.dr + J bVdn (2.11) 

r a r n 


That is, the work done by the forces of the first system with the displacements of 
the second system is equal to the work performed by the forces of the second 
system on the displacements produced by the first system. Now let the first system 
be the one we are seeking a solution to with the assumption that the body forces are 
neglected and let the second system correspond to the fundamental solutions for the 
traction and displacement due to a unit concentrated load fj in an infinite body, as 
shown in Fig.2.2. That is. 


b. = 0 

j 

t- = t.(Q) 
j j ^ 

Uj = u.(Q) 
b* = fj = 5(P,Q) e. 
tj =T..(P ,Q)e. 

u * = U..(P,Q) e. ( 2.i2) 
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Figure 2.2 Generalized Force Systems: (a) First System with Traction tj and 
Displacement u • (b) Second System of Kelvin’s Solution 
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Substitution of Eq.2.12 into Eq.2. 1 1, leads to the following equation 


Ju ij (P,Q)t.(Q)e.dr(Q) =jT ij (P 1 Q)a.(Q)e.dr(Q) +JS(P,Q)e j u j (Q)dQ(Q) 
r f a (2.13) 


Noting that 


J SCP^ejU.C^d^CQ) = u.(P) e. = u j (P) 5 . j e. 
n 


(2.14) 


and using Eq.2.14, Eq.2.13 can be rewritten as 


J U ij (P 1 Q)t j (Q)dT(Q)c i =jT ij (P,Q)u.(CBdr(Q)c. + (2. 15) 

r r 


or 
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(2.16) 


S^P) =J U ij (P ) Q)t j (Q)dT(Q) -J T.j(P,Q)Uj(Q)dT(Q) 

r r 


So far the point P is inside the boundary T. In order to have the equation relate 
only the points on the boundary surface F, we need to move the point P toward the 
surface. However, due to the singular nature of the kernels Ujj(P,Q) and T\j(P,Q) 
as r tends to zero, a limiting process must be employed in order to obtain the 
boundary integral equation. Let us first choose a new boundary 


r n = r_+r e 


where r £ is a surface of semi-spheric shape and T. is the rest of the surface as 
shown in Fig.2.3. F £ should envelop the point P such that P is still in domain Cl 
and thus Eq.2.16 is still available. With the new boundary, Eq.2.16 now becomes 


S..u.(P) =Ju..(P,Q)t.((^dr(Q) +Ju ij (P,Q)t.(Q)dT(Q) 

r - r e 

-jT ij (P,Q)u.(Q)dT(Q) -jT ij (P,Q)u j (Q)dT(Q) (2. 17) 

r- r e 
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As 8 tends to zero, the boundary T. becomes the original boundary T. We also note 
that when point Q is within the region f g we have 


dT(Q) = 8^sin0 d9 d0 
r = e 


r.i 



= ni 

f sinQ cos <|> 

: sin9 sin <f> 
l cos9 


Vm - 1 


r,i nj - r,j n* = 0 

tj(Q ) =tj(P) and uj(Q) = uj(P) 


(2.18) 


where tj(P) and uj(P) are the tractions and displacements at point P which are 
constants over the surface of the sphere as s tends to be zero. Substituting these 
relations into Eq.2.17, the second term on the right hand side of Eq.2.17 becomes 


0 9 

JUjj(P.Q)tj(Q)dr(Q) =JJ___a___[(3-4v)8^+n i n j ] ^sinOdSd^tCP) 

"e 0 0 


= _£Ljj[(3-tv)S..«.n^ in9ded(|) ( , 19) 

0 0 
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Figure 2.3 Domain and Boundary Used for Limiting Process 
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As e -» 0, the limit of the second integral on the right hand side of Eq. 2.17 
becomes 


UmJu ij (P,Q)t.(Q)dr(Q)=0 


( 2 . 20 ) 


The integral of Eq.2.19 is no longer singular and thus Eq.2.20 tends to be zero as 8 
is approaching zero. As for the fourth term on the right hand side of Eq. 2.17 


4> 0 


[t (P,Q)u (Q)dT(Q) = f f — — — [(l-2v)5.. + 3n n ]e 2 sin0d0d<!m.(P) 

J « J JJ sxd-vie 2 1J 1 J J 


0 0 


=8^S) JJI (1 ' 2v)5 ii + 3n - n il sined9d4 ’ (2 - 21) 


$ 0 


0 0 


Substimting Eq.2.18 into Eq.2.21 leads to a matrix expression for the kernel inside 
the integration and taking the limit as e approaching zero 
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$ 9 

lim f f T..(P,Q)u.(Q)ds(Q) 
£ -> o j j y j 

0 0 


-U.(P) 

87t(l-v) 


<t> 9 

(l-2v)+3sin 2 9cos 2 <j>, sin 2 0sin<{)cos<J>, sin9cos0cos<t> 

\(?) 


{j 

sin 2 9sin(j>cos4), (l-2v)+3sin 2 0sin 2 <j), sin0cos0 sin<j> 

u 2 (P) 

sin0d0d0 

0 0 

sin0cos9cos<j), sin0cos9sin<j>, (l-2v) + 3cos 2 9 

UjO*)' 



( 2 . 22 ) 


The integral ranges will depend on the local geometry of the point P. If point P lies 
on a smooth surface. 


0 = y and <j> = 2 tc 


and each term inside the matrix in Eq.2.22 becomes 




-1 
2 ’ 

0, 

•1 

"^(P)" 

lim fT..(P,Q)u.(Q)dr(Q) = 
£—►01 y j 


0 , 

-1 
2 ’ 

o 

1 

u 2 (P) 

j 

Pe 


0, 

0, 

_u 3 (P)_ 


(2.23) 
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or in Cartesian tensor form 


J L m o J T 8 (P.Q)u j (Q )dTCQ) = 4 5.. u.(P) 
r £ 


(2.24) 


Substituting Eq.2.20 and Eq.2.24 into Eq.2.17 leads to the boundary integral 
equation for P lying on a smooth surface 


I 5 ij u / p ) = / U i/P-Q)t J (Q)dr(Q) -jT ij (P J Q)u j (Q)dT(Q) 
r r 


(2.25) 


If P is not on a smooth surface, the boundary integral equation is 


C.. 

y 



TyfP.^u.CQ) dT(Q) 


(2.26) 


where 
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(2.27) 


Cy(P) Uj(P) = +l e imj Tj.(P,Q)dr(Q) 


which is only a function of the local geometry in the vicinity of point P. The 
integration in Eq.2.27 can be carried out by using Eq.2.22 with the appropriate 
integral range according to the local configuration of point P. The Cij(P,Q) can also 
be calculated from the concept of rigid body motion[24]. When the body 
undergose a rigid body translation, the surface is free of traction and the 
displacement is an arbitrary constant. By setting 


tjCQ) = 0 

uj(Q) = constant 


Eq.2.26 becomes 



(2.28) 


Comparing Eq.2.27 and Eq.2.28, it is apparent that the former depends on the local 
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geometry of each point and is different point to point. Thus it is very tedious to 
calculate. On the other hand, the latter is an integration on the whole boundary 
surface with a different kernel function and is easy to calculate for different points 
on the boundary. In fact, Eq.2.28 is merely a by-product of the second term on the 
right hand side of Eq.2.26 when it is carried out by numerical integration method. 
Therefore, with little effort, the term Cjj(P) can be easily obtained. Eq.2.28 is thus 
adopted in this paper to calculate the Qj(P) term. 

Eq.2.26 is also known as Somigliana's identity for three dimensional linear 
elastostatics with zero body force. In a well-posed boundary value problem, either 
traction or displacement in a direction on a boundary will be prescribed. Therefore, 
any corresponding unknown value on the boundary can be solved by the boun dar y 
integral equation of Eq.2.26. 

After the unknown boundary tractions and/or displacements have been 
solved, the displacement of any point inside the body can be solved by Eq.2.16 
from the boundary data. In order to obtain the stresses for the interior points, 
Eq.2.16 is differentiated and substituted into Eq.2.7. This results in 


V P) = J D i jk C p -Q)t i (Q)<ir(Q ) -Js ijk (P,Q)u.(Q)dr(Q) 

r r 


(2.29) 


22 


where 


VP.Q) - *^jk U im'm^ P ’Q) + t‘[U ij . k (P.Q) + U^/P.Q)] 


-1 

87c(l-v)r 


a-2v)( r , j 8 k i+ r, k 5 ir r,.8 jk ) + 3 r,s,x, k 1 (2.30) 


and 


S iik (P,Q) = XS^CP.Q) + utT^CP.Q) + T^oa 

=— Xj-j ( 3o m r, m [(l-2v)r. i 5 jk+ v(r. j 8 ]d+ r lt S ij )- 5 r.^r,,] 
+3v(n.r^+n k r, j )r,. + (l-2v)(3n.r,x, k + n.5 ki +n k 8. j ) 
- (l-4v)n.S jk } (2.31) 
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2.3 Numerical Implementation of Boundary Integral Equation 


Owing to the difficulties of solving the boundary integral equation in closed 
form for practical problems with complicated geometries and boundary conditions, 
it is necessary to solve the boundary integral equation numerically. 

The procedures for obtaining a numerical solution of the three dimensional 
boundary integral equation starts out with a discretization of the boundary surface 
T into m piecewise isoparametric surface elements. Such elements have been well 
developed for the Finite Element Method [25], Each element consists of n m nodes, 
the number depending on what kind of interpolation is employed in each element. 
The elements implemented in the present work include three-node linear triangular 
elements, four-node linear quadrilateral elements, six-nOde quadratic triangular 
elements, and eight-node quadratic quadrilateral elements. The shape functions for 
the isoparametric elements are derived in Appendix A. The cartesian coordinate of 
each node is given by 



(2.32) 


where xj a is the i - Cartesian coordinate of node a, N 0 ^) is the shape function for 
node a which is a polynomial function of intrinsic coordinates £ = ( ^ )• 

Eq.2.32 represents a one-to-one mapping of any point on the element from the 
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three dimensional Catesian coordinate space into the two dimensional ( q^ ^ 2 ) 
coordinate system, as shown in Fig.2.4. The isoparametric elements also use 
identical shape function to interpolate any function on the element, such as the 
displacements and tractions in the Boundary Element Method. Symbolically, 



(2.33) 


where Uj a , q a are the nodal values of displacement and tractions, respectively. By 
discretizing the surface T into m segments and utilizing the shape function, 
Eq.2.26 becomes 


c u “f Till Ju il (P.Q©)N a ©J 1 (5)drffif 

r, 

m ^ni C rt rrl 

- i Z ( E i jT ij (P > Q(^))N a ©J 1 ©dr©u al , P not summed (2.34) 
1*1 


where Cjj p and uj p stand for the coefficients Cjj and displacement in j-direction of 
nodal point P, respectively. The integral over the whole surface is carried out by 
summing up the integral over each element surface Tj . The tj 0 ^ and uj 0 ^ is the 
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Figure 2.4 Mapping of Isoparametric Elements 
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traction and displacement in i-direcdon of node a in element 1. The term Jj(^) is 
the determinant of the Jacobian matrix of element 1 deduced from the change of the 
integration variables from the Cartesian coordinate to the intrinsic ^-coordinate 
system. It is found to be 



3x. 


' sq 


a=l 0^ 1 a=l d $ 2 '1 

Ifi) = 



3x 2 9x 2 
’ *1 


(^ r ^ 2 ) ^ al 9N (^ r ^ 2 ) ^ al 

a=l 0^ '2 ’a=l 3^ 2 


(2.35) 

Similary, Eq.2.28 is rewritten as 


^ij l T..(P,Q(5))J 1 (5)dT(5) 
l*i 


(2.36) 


A Gaussian quadrature scheme can now be applied to accomplish the integration. 
For the case when the point P is not on the element which is being integrated over, 
the Gaussian quadrature is straight forward, that is. 
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J>F,(p,ydr© =| ib | w; w Jp l( p, 5 | 4 j) (2.37) 


where 


'P 1 (P^)=U ij (P,Q(^))N a (^)J 1 (^) 


or 

=T ij (P,Q(^))N a (^)J 1 (^) 


In Eq.2.37 n a and n^ are the order of Gaussian quadrature and Wj a and w ? b are 
the weights of the corresponding Gaussian integration points ^ ^ > respectively. 
When point P is one of the nodes of the element under integration, the stan dar d 
Gaussian quadrature will not give accurate results because of the 1/r and 1/r— 
singularities of Uy (P,Q) and Tjj(P,Q), respectively. Therefore, special treatment of 
this singular integral must be used in order to obtain an accurate solution. The 
method employed here follows the work of Rizzo and Shippy [15]. The element is 
divided into triangles according to the position of the point P as shown in 
Appendix B. Inside each triangle, the ^-coordinate system is transformed into a 
local polar coordinate system , r and 9, and the term in the integral becomes 
rdrdG. The additional r due to this transformation eliminates the 1/r singularity of 
the kernel Uy(P,Q). Gaussian quadrature can then be applied to the polar 
coordinate system. Of course, one more transformation is needed to map the r,0 
coordinate system to another polar coordinate system r and £ so that the range 


28 



varies from -1 to 1. The details of the procedures are illustrated in Appendix B. 
For the integration of the kernel Ty(P,Q), a 1/r singularity remains after the 
transformation from the ^-coordinate to the polar coordinate. However, by 
substituting Eq.2.36 into Eq.2.34, Eq.2.34 can be rewritten as 


5 jT ij (P,Q(^))J 1 (q)dr(^)u p 
r i 

=| lc | 1 Ju ij (P,Q(^))N a (^)t“ l J 1 (^)dr(^) 

r i 


It Jt J 

Pi 

(2.38) 


Distinguish the elements containing point P from those elements without point P for 
the integration involving the kernel Ty(P,Q), and Eq.2.38 becomes 


3* jT. j (P > Q(^))J 1 (^)dr(^] > 

Pi 

=33 Ju ij (P,Q(^))N a (^)t“ 1 J 1 (^)dr©- 33; jT. j (P,Q(^))N a (g)u“ 1 J 1 (^)dr(^) 

Pi r l 

m p n m r n ^ m p r p 

-2^ j T..(P,Q(^))N (^)u“ 1 J 1 (^)dT(^)+ 1 I jT. j (P,Q(^))J 1 ©dr(q)u[ (2.39) 

Pi Pi 
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where m x is the total number of the elements excluding point P and mp is the total 
number of the elements including point P. Obviously, m = m x + mp. Combining 
the last two terms on the right hand side of Eq.2.39, results in 


III J U ] /P,Q©)N“©J 1 ©dr©t'‘ 1 jT (j (P,Qffl)N“® ^©dT ®^ 1 

r t r. 

- ,jl, Jr ij (P>Q©)<N a © -s^)J 1 ®dr©u» 1 (2.40) 

r, 


where 


S = J T yOP.Q®)dr© 


(2.41) 


Eq.2.41 is the negative of the integral of the kernel Tij(P,Q) over the elements 
which do not contain the point P, hence, Eq.2.41 is no longer singular. As for 
Eq.2.40, the remaining 1/r singularity of the kernel TjjCP.Q) is removed by the 

special shape function N 0 ^) - 8 a p in the last term on the right hand side of 
equation because when a coincides with point P ,the constant term is eliminated 
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and thus the smallest order of the shape function is of order r. This additional r can 
cancel the 1/r singularity after the coordinate has been transformed into the local 
polar coordinate. Using this procedures all the singularities are eliminated and 
Gaussian quadrature can be used. 

Eq.2.40 represents an equation at each discrete point P on the boundary 
surface constraining the boundary displacements and boundary tractions in the i- 
direction. For a surface including N nodes, Eq.2.40 can be expressed as 


J 


— p P 
C.u. - 
y j 


u N 

V. J J 


r ^ 
J 


'12 N 


U J 


1 2 

N 1 

T.,T., X- 

< 


> = 

U...U.,... 


_ — y — u *-yj 

1 

1 • 1 

[" 

— 1 J — 1J 



< 


2 

t. 

J 


r (2.42) 


; n 

V. J J 


where each term of Xij for node n is the integration summed up from the 

contribution of the elements which share the same node n. The same applies for 
n 

. Eq.2.42 can be rewritten as 


C P u p + 
y j 




(2-43) 
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This equation can be expressed more simply by combining the first term into the 
summation on the left hand side, i.e.. 


N *n r> 1^ *n n 

I T. u = I U. t 

n=l y J n=l >J J 


N 


(2.44) 


where 


T* n = £ for n * P 

ij 

T ?“£i + C ij forn = p (2.45) 


Another way of dealing with the singularity which is worthy to mention here is 
through the use of rigid body motion. For a rigid body motion in direction j, 
Eq.2.44 reduces to 


N « n 

I T =0 
n=l y 


(2.46) 


This indicates that the sum of Tjj* n from each node in a row for certain direction j 
should be zero. Hence the value of the singular term Ty*P when n = P can be 
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easily calculated once all the other terms are known, that is. 


*p N «P 

T. =- 1 T. (2.47) 

U n=l ‘J 

n*P 


For three dimensional problems the indices i and j range from 1 to 3. Therefore, 
for a surface including N nodes, the dimension of the total algebraic system of 
equations formed by Eq.2.44 for each node in each direction is 3Nx3N. The 
system can be represented in matrix form as 


T*u = U*t 


(2.48) 


The ma trixes T* and U* are rearranged by interchanging the suitable columns on 
each side of equation so that all the unknown variables are contained in a column 
vector x and all the prescribed values of the boundary are included in the column y 
on the other side of the equation. Symbolically, 


Ax = By = f 


(2.49) 


The system of equations can now be easily solved by the Gaussian elimination 
method. 
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After all the boundary data are obtained, the displacements and stresses at 
interior points are solved by substituting the boundary tractions and displacements 
into Eq.2.16 and Eq.2.29, respectively, and using the shape functions to carry out 
the integration. 


2.4 Subregion Technique 


Since the boundary integral equation is a constraint relation between 
boundary tractions and boundary displacements with the kernel functions which 
include the term 1/r and l/r^, any two distinct points on the boundary can not 
coincide. Therefore problems for which the boundary includes two contacting 
surfaces can not be solved by the Boundary Element Method using a single region. 
Partitioning the whole boundary into subregions is necessary to deal with such 
problems. 

Consider for simplicity the case where the body is partitioned into two 
subregions. Similar procedures can be followed to separate the body into more 
than two regions. A body with domain Q surrounded by the boundary T is 
partitioned into two subregions; one consists of domain Qj and boundary Tj and 
the other possesses the domain Q 2 an d boundary ^ , as was shown in Fig.2.5. 
The two regions share the same interface Tj. Each subregion can be treated as an 
independent body. Thus the procedure described in Section 2. 1 can be used to 
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Figure 2.5 Domain Divided into Two Subregions 
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form the system of algebraic equations for each subregion. Before the compatibility 
and equilibrium conditions are enforced along the interface between two 
subregions, the corresponding columns of the matrix U* of unknown tractions and 
the corresponding columns of the matrix T* of known displacements in Eq.2.48 
must be interchanged so that all the known and unknown variables will be at the 
same side of the equations. Let A s and B s be the kernel matrices related to the 
unknown column vectors Uk g s and K e * on the external surface of region s after 
rearrangement. Furthermore let tj s , uj* be the unknown tractions and unknown 
displacements on the interface of region s. The system of equations for subregion 
one can be expressed as 




Ukl 



A 1 


L 1 

1 


B 1 


(2.50) 


After pre-mul tiplying each side by the inverse of Bj, Eq.2.50 becomes 


dL 

d L 

d L 

D 1 . 

U 




L 1 

1 


£ 


ur 


(2.51) 
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where A[. By the same procedure the system for the second subregion 

becomes 


7 

Dec 

D 2 . 

ei 


D 2 

U 


uk: 


t 2 

L 


k: 


u 2 

l 


(2.52) 


Enforcing the compatibility and equilibrium condition along the interface 


ti 1 =-ti 2 = ti 


1 4 

ui i = ui =ui 


(2.53) 


Equations 2.51 and 2.52 can be combined to be [26] 


dL 

d L 

0 

d L 

1 

a 

D 2 

ie 

0 

-D 2 . 

ei 

d 2 « 


Ul4 


Uk! 


Ki 


0 


K* 


(2.54) 
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(2.54) 

This is a system of 3( N e ^+ N g 2 + Nj ) equations for 3N e ^ unknowns on the 
external surface of subregion one, 3N g 2 unknowns on the external boundary of 
subsurface two, and 3Nj unknown tractions on the interface. The unknown 
displacements on the interface can be obtained through either Eq.2.51 or Eq.2.52 
after Eq.2.54 has been solved. The equation which is not used can be used as a 
check. 


This method is simple and direct. However it is not usually adopted 
because the inverse of the matrix of kernel B s must be calculated for both 
subregions. This involves a tremendous amount of computer processing time and 
requires a lot of memory space. Therefore, this method is not used in this paper. 

An alternative procedure which eliminates the need for solving the inverse 
of the matrix is described herein. After the appropriate columns of T* and U* 
have been interchanged, the system of equation of the two subregions are 
assembled into two big matrices and two col umn vectors as 


A 1 

0 

0 

*7 

A" 



1 


Uk e 2 


1 


B 1 

0 

0 

B 2 


K 


u. 


k: 


u 2 


(2.55) 
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To be more specific, 



bL 

B 1 . 

ei 

0 

B 1 

ie 

B 1 . 

U 

0 

Bee 

B 2 

ei 

B-c 

B 2 

a 




1 


(2.56) 


Using the compatibility and equilibrium conditions (Eq.2.53), Eq.2.56 reduces to 



(2.57) 


Since now the known variables appear on the right hand side of the equation and 
the unknown variables appear on the other side,the equations can be solved without 
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any difficulty. Although this method avoids solving for the inverse matrix and thus 
saves a lot of computer processing time, it still has the drawback of using two big 

matrices which occupy tremendous storage space. Therefore another algorithm is 
introduced next. 

Let the numbering of the nodes along the interface be always arranged at the 
position following the external nodes for both subregions, as shown in Eq.2.50. 
Rewrite Eq.2.50 for the first region as follows 


A 1 

a! 

e 

i 


UK 


n) 

l 


B 1 

B 1 

e 

1 


K 


(2.58) 


where the dimensions of each block axe as following 


A* and B* are 3(N* + N.)x 3N* 
Aj and B* are 3(N* + N.)x 3N. 

Kg and Uk] are 3NgX 1 
t. 1 and u. 1 are 3N : x 1 

l 1 l 


Bg 1 and IQ 1 can be multiplied together to form a known column vector f 1 
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Al 

A 1 

l 


Ukg 


f 1 

+ 

B 1 

1 




u 1 

1 






(2.59) 


Moving the second term of the right hand side of Eq.2.58 to the left hand side leads 
to 




(2.60) 


Note that the left hand side of Eq.2.60 now become a matrix with dimension of 
3(N e l + Nj) by 3(N e l + 2N0 multiplied by an unknown vector with dimension of 
3(N e ^ + Ni) by 1. Obviously, this system can not be solved since there are 3Nj 
more unknown than the total number of equations. Before we proceed to the 
second subregion, Eq.2.60 is reduced to 



(2.61) 


41 





Eq.2.61 can be represented more conveniendy as 



(2.62) 


Eq.2.62 can be treated as two set of equations ,viz. 


X 



+ 

Gl 


u. 1 

1 

t. 1 

zz 

g 1 





1 




(2.63) 


Hi 




(2.64) 


The same procedure can be applied to the second subregion and a similar set of 
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equations is obtained 




i 


( 2 . 66 ) 


Note that the unknowns in Eq.2.64 and Eq.2.66 now only involve the unknown 
displacements and unknown tractions on the interface between the two subregions. 
By applying the compatibility and equilibrium condition of Eq.2.53, 


t i 1 =-t i 2 = t i 

1 2 
u r = u i 


these two sets of equations can be assembled together as Eq.2.67 to solve for the 
unknown tractions and unknown displacements along the interface. 
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Hi 


u. 

i 



H 2 


t. 

1 


h2 


(2.67) 


where 



After the displacement uj and the traction tj of the interface have been solved using 
Eq.2.64, the unknowns on the external surfaces of each subregion can be solved 
through Eq.2.63 and Eq.2.65. 

This method not only saves computer excution time by eliminating the need 
for solving for the inverse of the matrices but also saves a lot of memory storage 
space because the procedure can be performed on the system of equations of each 
subregion separately with little effort by using an index control in the computer 
program. It is therefore used in this paper. 
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2.5 Discontinuity of Traction 


The data of prescribed tractions and prescribed displacements must be read 
into the computer according to the element rather than the node due to the possibility 
of a node lying on a position without an unique tangent plane. That is, when a 
node is at the intersection of two or more planes, the traction acting on this node 
depends on which plane is considered, since each surface is associated with a 
different normal. For example, in Fig.2.6, three elements share the same node P. 
For the sake of easy interpretation, we let the unit normal of these three elements 
coincide with the unit vectors of the Cartesian coordinate system. 


n 


1 


A 

= 1 


= J 



The component of the traction of node P in direction i on element 1 is the normal 
stress a} i but the traction of point P in the same direction on element 2 is the shear 
stress aj 2 me! the traction in the i-direction on element 3 is the shear stress at 
point P. Therefore, the prescribed traction must be input according to the element 
so that each different traction on the same node can be multiplied by an appropriate 
kernel contributed from each element surface. Also note that when the displacement 
in one direction is prescribed at a certain node, for a unique solution only one of 
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the corresponding tractions on this node can be unknown and all the others must be 
prescribed. For the same example in Fig.2.6, two of three tractions on three 
different elements must be known if the displacement of that node is prescribed. 
Eq.2.37 in fact is a concise expression for convenient intrepretation. To be more 
specific, it should read 



? 


? 

? 


(2.68) 


Thus, the equation for a certain row is 


tV 


■uft + u £ + u,V + 


(2.69) 


When the tractions tqP, t 2 P and t 3 p are known at point P, they must be multiplied 
by these different kernel function U^, and U 3 P which are calculated from 
each different element with different normal vector. As for the case of the 
displacement being prescribed, only one traction can be unknown, say t 2 ^, and the 
traction on another surface t]P and t 3 ^ must be forced to be known. Thus, the 
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Figure 2.6 Illustration of Different Tractions at Node P 
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equation becomes 



so that only one unknown is left on the right hand side of the equation which is 
needed to be solved and all the known values are Summed up on the left hand side 
of the equation. 

One incident merit of inputting the traction element by element is that 
discontinuous tractions can be modeled exactly. 
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CHAPTER THREE 


LINEAR ELASTIC FRACTURE MECHANICS 


3.1 Introduction 

The presence of cracks in a material causes a stress concentration in the 
vicinity of the crack tip. Consequently, plastic yielding or microcradting will occur 
in the region surrounding the crack tip. Linear elastic fracture mechanics assumes 
that the nonlinear deformations are restricted to a region whose dimensions are 
small compared to other characteristic dimensions, so that the elastic solution 
provides an accurate description of the stress and displacement fields in the vicinity 
of the crack tip. This is often referred to as small scale yielding (S.S.Y.). 

It can be shown that the loading on a crack is in general a superposition of 
three independent modes as shown in Figure 3.1. The first (Fig.3.1.a) is called 
the crack opening mode, or mode-I, which is a result of a relative normal separation 
of the crack surfaces (symmetric with respect to x-z and x-y planes). Fig.3.1.b is 
called the crack sliding mode, or mode-II, which is associated with a relative sliding 
displacement in the x-direction (symmetric with respect to x-y plane and skew- 
symmetric with respect to the x-z plane). The tearing mode (mode-EH) corresponds 
to relative motion in the z-direction of the two crack surfaces (skew-symmetric with 
respect to the x-y and x-z planes). Using Westergard's technique, Irwin and 
Williams showed that the stress and displacement fields in the vicinity of a crack tip 
can be expressed as an infinite series whose leading term is square root singular. 
The coefficient of this singular term is defined as the stress intensity factor. With 
respect to a r-0 polar coordinate system as shown in Figure 3.2 , the stresses and 
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Figure 3.1 The Three Basic Modes of Crack Surface Displacement 
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displacements near a crack dp are given by [27] 


Mode I 


K. 


hr cos|- [ 1- sin-f- sin^r ] 


^xx , 

V27CT - ^ 

<Jyy= —= cosy [ 1+ siny siny- ] 
V27CT 2 22 


K 


a xy~ 


1 sinyCOSyCOS-y 


V 27tr 

^zz = ^ (^xx"*" ^yy )* ^xz~ ^ yz~ ® 
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U 


U 


1 


2 


3 



cosy [ 1 - 2v + sin 2 y ] 
sin-y [ 2 - 2v - cos 2 y ] 


(3.1) 


Mode II 


a xx =y==T siny [ 2 + COSy COSy- ] 

■Jim 2 2 2 

a yy = siny COSy COSy 

V 2tct - z - 
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K. 


u i = 


n 


P- 

K n 
U 2 u 


V^ 5i "f [2-2 v + C0 S 2 |j 


cos-f- [- 1 + 2v + sin 2 -^- ' 
2tc 2 2 ' 


u 3 = 0 


(3.2) 


Mode EH 


0x2 = ' “'z 


a xx _ ®yy a zz~ ^xy ~ 0 



(3.3) 


where Kj, Kji, and Kjjj are the stress intensity factors corresponding to modes I, 
13, and HI, respectively. Note that Eq.3. 1-3.3 are valid only when r « L, where L 
is another characterictic length of the geometry (may be the crack length). Also note 
that these equations do not contain any information about the externally applied 
loading, crack geometry, geometric configuration, etc, and that the stress intensity 
factors are not functions of the local coordinate r and 8. These factors are embedded 
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Figure 3.2 Coordinate and Stress Components for Crack Tip Stress Field 
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in K[, Kq, and Kjjj That is, the stress intensity factors Kj, Kjj, and Kjjj, which 
must be determined by the boundary conditions and the loadings control the 
magnitudes of the stress and displacement fields around the crack tip. Thus, a 
given combination of values of Kj, Kq, and Kjji represents uniquely a crack tip 
stress field environment for small scale yielding. The determination of the stress 
intensity factors is thus the most important task in linear elastic fracture mechanics. 
Although Eq.3. 1-3.3 are valid for the plane strain problem, they may be modified 
to represent the plane stress problem by letting a z = 0 and substituting v with 
v/(l+v). 

Eq.3.1 and Eq.3. 2 were derived for plane problems. However, it was 
shown by Sih and Liebowitz [28] that for an elliptical crack in a three dimensional 
linear elastic body, the local stress and displacement fields along a crack front are a 
superposition of plane strain and antiplane shear. Hence, Eq.3. 1-3.3 can still be 
used for three dim ensional crack problems as long as the coordinate system is 
allowed to move along the crack front with its z-axis tangent to the crack front and 
the y-coordinate perpendicular to the crack surface, as shown in Fig.3.3. Referring 
to the moving coordinate, the stress intensity factors are now also a function of the 
position of the origin on the crack front and the formulas are not available for a 
crack near a free surface because the stress singularity is not of inverse square root 
there. 
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Figure 3.3 Coordinate System on the Crack Front for Three Dimensional Crack 
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3.2 Calculation of Stress Intensity Factors Using the Boundary Element Method 


The variation of any function in an isoparametric element is polynominal. 
The degree of the polynominal depends on the number and arrangement of the 
nodes. Therefore, if a quadratic surface element is used in the vicinity of a crack, 
the distributions of the displacement and traction in the element will have at most 
quadratic variations. Since the variation of the displacement is square root of r and 
the variation of the tracton is inverse of square root of r around the crack, a fine 
mesh is needed to model the crack so that the quadratic variation can imitate the 
correct distribution inside each small segment. However, even this refinement 
cannot achieve a high degree of accuracy. 

Fortunately, this problem was solved by Barsoum [18] who modified the 
quadratic isoparametric element by relocating appropriate midside nodes to the 
quarter-point to capture the inverse square root singularity. Even though he had 
done this for the Finite Element Method, a similar approach can be used in the 
Boundary Element Method. For example, consider an eight node quadrilateral 
element with two sides having equal length L perpendicular to the crack front ( side 
1-5-2 ),as shown in Fig.3.4. Relocate the two midside node ( 6 and 8 ) to the 
quarter-point near the crack front. Denoting the distance originated from the crack 
front to any point on the element by r, then 


N i^r^ r i 


(3.4) 
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Figure 3.4 Illustration of 8-node Quarter Point Element 
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Choosing ^ - 0, = L, and = rg = L / 4, Eq.3.4 becomes 


r - \ (5, + V 1X1+ ?,)(!+ 5j) L + i (5,- 1X$,- l)(5i +1) L 

+i (1+ 5,X1 - £)xi- + 1 (1+ yd - $L +i(l - ?,xi • 


Simplifying, 



(3.5) 


Substitution of Eq.3.5 into Eq.2.33, leads to the following variation of the 
displacement and traction in direction i versus the distance r 
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u.(r) = l N k u k 

i k=l i 

= [(-u! - u 2 - u? - u 4 - u 5 + 2u? + u 7 + 2u 8 ) 

l l 1 1 1 1 1 l ' 

+ 4 1 (3uf - 3u 2 - u 3 + u? + 4u? - 4u? ) 

2 i ii i i i 

, 2 

+ — - (- u* - u? + u? + u 4 + 2u 3 - 2u 7 )] Af- 
2 1 1 1 1 1 1 v L 

+ [( u. 1 + u 2 + u 3 + u 4 - 2u 6 - 2u 8 ) 

1 l 1 l 1 1 ' 

E E 2 

+ 1 <4 * y<- u l + u r > + t< u ! + “? • 2u f > i p 6) 


which can be rewritten as. 


u.(r) = a! + A 2 /f + A 3 f 
i i iV l 1 L 


(3.7) 


Similary, 


t.(r) = B 1 + B 2 U- + B 3 f 
i i i V L i L 


(3.8) 
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The displacement variation given in Eq.3.7 contains the square root of r term which 
is asymptotically correct. However, the traction appearing in Eq.3.8 does not 
include the correct inverse square root of r term needed to model the stress 
singularity around the crack tip. The correct singularity is obtained by multiplying 
the right hand side of Eq.3.8 by the factor [19] 



such that the variation of traction becomes 

•iW = ( b[ + + B i 

- b !V? +b ? +b iV? (3 - 10) 


which possesses the correct inverse square root of r term. 

For six node triangular quadratic elements, if all the sides are straight as 
shown in Fig.3.5, by similar procedure, the relation between r and the nature 
coordinate ^ can be obtained as 
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Figure 3.5 Illustration of 6-node Quarter Point Element 
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(3.11) 


where c is the ratio 1 S 3 which is constant for a given direction r. It is easy to 
show that the variation of the displacement has the same form as Eq.3.7 by 
substituting Eq.3.11 into Eq.2.30 with the help of appropriate shape function. The 
correct traction variation is obtained by multiplying the shape function for traction 
by the correction factor 



(3.12) 


The stress intensity factors are then evaluated by the displacement correlation 
technique [29], By setting 0 = 180° in Eq.3.1 -3.3, the displacement fields become 


u i = 


U 2 = 
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(3.13) 
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Equating the square root of r term in Eq.3.13 to that in Eq.3.6 yields 


K r = 


2(1 -v 2 ) 


^7 [(-Uj 1 - Uj 2 - Ul 3 - Uj 4 - Ul 5 + 2 Ul 6 + u t 7 + 2u x 8 ) 

+ j 1 (3UJ 1 - 3 u t 2 - + u x 4 + 4u : 6 - 4u x 8 ) 

? 2 

+ Y (- u x ! - u 2 + Uj 3 + n 4 + 2Uj 5 - 2UJ 7 )] 


K n " 


2(1 -v 2 ) 


t [ ^2 * U 2 2 - U 2 3 * U 2 4 * U 2 5 + 2U 2 6 + U 2 7 + 2U 2 8 ) 

+ fou , 1 - 3u 2 2 - + u 2 4 + 4u 2 6 - 4u 2 8 ) 

+ -± (- U 2 1 - U 2 2 + u 2 3 + U 2 4 + 2u 2 5 - 2U, 7 )] 




2(l+v) 


[(-U3 1 -u 3 2 -U3 3 -U3 4 -U3 5 + 2u 3 6 + U3 7 + 2u 3 8 ) 

+ |- 1 (3u 3 1 - 1^ 3 + i^ 4 + 4u 3 6 - 4^ 8 ) 

E 2 

+ ^-(-u 3 1 -u 3 2 + u 3 2 + u 3 4 + 2u 3 5 -2u 3 7 )] 


(3.14) 
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For non-symmetric problems both crack surfaces must be modeled. The stress 
intensity factors are given by 




4(l-v) 


£ ([(- « 3 - u, 4 + 2 u, 6 + 2u, ! ) - (- u', 3 - u‘+ 2u> u* 7 + 2^)] 

+ 1[( - Ul 3 + Uj 4 + 4 Ul 6 - 4 Ul 8 ) - ( - u ; 3 + u 4 + 4 U ; 6 - 4 U ; 8 > ] 
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K m = 


4(1+v V % 3 - 2 V + 2 S 8 >-(- 2 ^ <?+ ^ 

+ -| l [(- U3 3 + U3 4 + 4U3 6 - 4U3 8 )-( - tlj 3 + u^ 4 + 4 i^ 6 - 4 u* 8 )] 


+ y [ ( U3 3 + U3 4 - 21X3 7 )-( u 3 * 3 + u^* 4 - 2U3* 7 )] } 


(3.15) 


where the asterisk refers to the displacement of the node on the element opposite to 
the one shown in Fig.3.5. Note that Eq.3.14 and Eq.3.15 involve the natural 
coordinate ^ up to quadratic terms. Therefore a quadratic variation of the stress 
intensity factor in the ^ direction can also be represented. 
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CHAPTER FOUR 


VERIFICATION PROBLEMS 

Several verification problems were studied to examine the accuracy and 
efficiency of the Boundary Element Method. This chapter presents the results 
obtained using the techniques described in the previous chapters for problems 
whose analytical solution is known as well as for a problem whose solution was 
obtained numerically by other researchers using different techniques. All problems 
were performed using single precision on a CRAY-X-MP. The material properities 
are: E = 30,000 ksi and v = 0.3. 

4.1 Prismatic Bar under Uniform Tension 

The first verification problem is a linear elastic prismatic bar under simple 
tension as shown in Figure 4.1. Only one-eighth of the specimen is modeled due to 
the symmetry of this problem. The domain is chosen to be a cube with side length 
equal to 1 subjected to a uniform tension on the surface z=l. The planes x=0, 
y=0,and z=0 are fixed in x, y, and z directions, respectively. The geometry and 
boundary segments are given in Figures 4.2 to 4.5. Four different elements: 3- 
node triangular element, 4- node quadrilateral element, 6-node triangular element, 
and 8-node quadrilateral element were tested. Different schemes of Gaussian 
quadrature were chosen to perform both the regular integral and singular integral as 
described in Section 2.3. The traction in direction z at point A and the transverse 
displacement in the x-direction of point B in Figure 4.1 were calculated and 
compared with the exact solution [31]: 
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Figure 4.1 Geometry and Boundary Conditions for Prismatic Bar in Uniform Tension 
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Figure 4.2 3-Node Triangular Element Mesh for Uniform Tension Problem 
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Figure 4.3 4-Node Quadrilateral Element Mesh for Uniform Tension Problem 
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Number of Nodes: 26 
Number of Elements: 12 


Figure 4.4 6-Node Triangular Element Mesh for Uniform Tension Problem 
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Number of Nodes: 20 
Number of Elements: 6 


Figure 4.5 8-Node Quadrilaterial Element Mesh for Uniform Tension Problem 
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t A = -a 


1^= -V 0/E 


(4.1) 


The results are summarized in Table 4.1 to 4.4. The results indicate that as long as 
the singular integral is evaluated with at least a 3x3 quadrature, then the accuracy 
depends on the number of points used to evaluate the regular integral. For reliable 
results the regular integral should be evaluated using at least a 3x3 quadrature for 
quadrilateral elements and 6 points quadrature for triangular elements. 

Tables 4.5 to 4.8 demonstrate the accuracy of the multi-domain technique. 
The conclusions are similar to those for the single region. However, the two- 
region results are more accurate since for the same number of integration points 
there are more nodes. 
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NRI 

NS I 

u A (><10' 6 ) 

Error (%) 

C B 

Error (%) 

CPU Time 
( second ) 

1 

lxl 


>100 


>100 

0.0628 


2x2 

1 .4397 

43.97 

-0.4897 

51.03 

0.1156 


3x3 

1 . 4214 

42.14 

-0.5409 

45.91 

0.2056 


4x4 

1.4197 

41 . 97 

-0.5451 

45.49 

0.3362 


5x5 

1.4197 

41.97 

-0.5454 

45.46 

0.4948 


6x6 

1.4197 

41.97 

-0.5454 

45.46 

0.7030 

3 

lxl 

0.7303 

26.97 

-0.4628 

53.72 

0.1110 


2x2 

0 . 9376 

6.24 

-0.9861 

1.39 

0.1629 


3x3 

0.9600 

4 . 00 

-1.0470 

4.70 

0.2430 


4x4 

0 . 9582 

4 .18 

-1.0611 

6.11 

0.3808 


5x5 

0 . 9580 

4.20 

-1.0611 

6.11 

0.5437 


6x6 

0.9580 

4.20 

-1.0611 

6.11 

0.7409 

4 

lxl 

0.7244 

27.56 

-0.4808 

51.92 

0.1312 


2x2 

0 . 9310 

6.90 

-1.0121 

1.21 

0.1852 


3x3 

0 . 9548 

4.52 

-1.0800 

8.00 

0.2762 


4x4 

0 . 9549 

4.51 

-1.0862 

8 . 62 

0.3998 


5x5 

0 . 9547 

4.53 

-1.0867 

8.67 

0.5635 


6x6 

0 .9547 

4.53 

-1.0867 

8.67 

0.7684 

6 

lxl 

0.7564 

24.36 

-0.4041 

59.59 

0.1954 


2x2 

0 . 9746 

2.54 

-0.9250 

7.50 

0.2466 


3x3 

1.0035 

0.35 

-0.9907 

0.93 

0.3379 


4x4 

1.0017 

0.17 

-0.9967 

0.33 

0.4635 


5x5 

1 . 0019 

0.19 

-0.9971 

0.29 

0.6254 


6x6 

1.0019 

0 . 19 

-0.9971 

0.29 

0.8280 

7 

lxl 

0 . 7559 

24.41 

-0.4031 

59.69 

0.2184 


2x2 

0 . 9739 

2.61 

-0.9237 

7.61 

0.2692 


3x3 

0 . 9984 

0.16 

-0.9897 

1.03 

0.3623 


4x4 

0 . 9988 

0 . 12 

-0.9956 

0.44 

0.4865 


5x5 

0 . 9989 

0.11 

-0 . 9960 

0.40 

0 . 6615 


6x6 

0.9989 

0 . 11 

-0.9960 

0.40 

0.8518 

12 

lxl 

0.7552 

24.48 

-0.4072 

59.28 

0.4131 


2x2 

0 . 9728 

2.72 

-0.9280 

7.20 

0.4764 


3x3 

0 . 9979 

0.21 

-0.9937 

0.63 

0.5624 


4x4 

1 . 0002 

0.02 

-0.9997 

0.03 

0.6848 


5x5 

1 . 0001 

0.01 

- 1.0000 

0.00 

0.8425 


6x6 

1.0001 

0.01 

- 1.0000 

0.00 

1.0503 


NRI : Number of Gaussian quadrature points for regular integral 
NS I : Number of Gaussian quadrature points for singular integral 
Exact solution : u=l.E-06, t=-1.00 


Table 4.1 Beam in Uniform Tension, 3-node^ Triangular Element: Comparison 
of Error for Displacement and Traction, and CPU Time 
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NRI NSI u a (*iq- 6 ) Error (%) t B Error (%) ^secMd? 

lxl 1 x 1 >100 >100 0.0565 

2x2 1.3438 34.38 -1.0000 0.00 0.1445 

3x3 1.3453 34.53 -1.0000 0.00 0.2951 

4x4 1.3453 34.53 -1.0000 0.00 0.4984 

5x5 1.3453 34.53 -1.0000 0.00 0.7735 

6 x 6 1.3453 34.53 -1.0000 0.00 1.0980 

2x2 lxl 0.9264 7.36 -1.0000 0.00 0.0764 

2x2 0.9559 4.41 -1.0000 0.00 0.1659 

3x3 0.9568 4.32 -1.0000 0.00 0.3120 

4x4 0.9568 4.32 -1.0000 0.00 0.5217 

5x5 0.9568 4.32 -1.0000 0.00 0.7866 

6 x 6 0.9568 4.32 -1.0000 0.00 1.1149 

3x3 lxl 0.9689 3.11 -1.0000 0.00 0.1118 

2 x 2 1.0008 0.08 - 1.0000 0.00 0.2003 

3x3 1.0017 0.17 -1.0000 0.00 0.3492 

4x4 1.0017 0.17 -1.0000 o.OO 0.5586 

5x5 1.0017 0.17 -1.0000 0.00 0.8286 

6 x 6 1.0017 0.17 -1.0000 0.00 1.1609 

4x4 lxl 0.9673 3.27 -1.0000 0.00 0.1582 

2 x 2 0.9992 0.08 -1.0000 0.00 0.2502 

3x3 1.0000 0.00 -1.0000 0.00 0.3997 

4x4 1.0000 0.00 -1.0000 0.00 0.6060 

5x5 1.0000 0.00 -1.0000 0.00 0.8679 

6 x 6 1.0000 0.00 - 1.0000 0.00 1.2027 

5x5 lxl 0.9673 3.27 -1.0000 0.00 0.2204 

2x2 0.9991 0.09 -1.0000 0.00 0.3112 

3x3 1.0000 0.00 -1.0000 0.00 0.4604 

4x4 1.0000 0.00 -1.0000 0.00 0.6677 

5x5 1.0000 0.00 -1.0000 0.00 0.9362 

6 x 6 1.0000 0.00 - 1.0000 0.00 1.2642 

6 x 6 lxl 0.9673 3.27 -1.0000 0.00 0.2985 

2x2 0.9991 0.09 -1.0000 0.00 0.3882 

3x3 1.0000 0.00 -1.0000 0.00 0.5335 

4x4 1.0000 0.00 -1.0000 0.00 0.7416 

5x5 1.0000 0.00 -1.0000 0.00 1.0112 

6 x 6 1.0000 0.00 - 1.0000 p . qq 1.3333 

NRI : Number of Gaussian quadrature points for regular integral 
NSI : Number of Gaussian quadrature points for singular integral 
Exact solution : u=l.E-06, t=- 1 .00 


Table 4.2 Beam in Uniform Tension, 4-node Quadrilateral Element: Comparison 
of Error for Displacement and Traction, and CPU Time 
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NRI 

NS I 

u^xiO ' 6 ) 

Error (%) 

C B 

Error (%) 

CPU Time 
( second ) 


1 

lxl 


>100 


>100 

0.3332 


2x2 

1 . 8007 

80.07 

-0.4114 

58.86 

0.5378 


3x3 

1 .8843 

88 . 34 

-0.3102 

68.98 

0 .8721 


4x4 

1 . 8754 

87.54 

-0.3153 

68.47 

1.3492 


5x5 

1.8759 

87.59 

-0.3239 

67.61 

1.9639 


6x6 

1 . 8776 

87.76 

-0.3224 

67.76 

2.7301 

3 

lxl 


>100 


>100 

0.5478 


2x2 

0.7482 

25.18 

-1.3980 

39.80 

0.7635 


3x3 

0.8437 

15 . 63 

-1.2393 

23.93 

1.0886 


4x4 

0 . 8473 

15.27 

-1.2196 

21 . 96 

1.5714 


5x5 

0 . 8459 

15 . 41 

-1.2299 

22.99 

2.1947 


6x6 

0.8465 

15.35 

-1.2291 

22.91 

2.9585 

4 

lxl 


>100 


>100 

0 . 6481 


2x2 

0.7151 

28.49 

-1.4042 

40.42 

0.3548 


3x3 

0.8083 

19.17 

-1.2411 

24.11 

1.1982 


4x4 

0.7938 

20 . 62 

-1.2162 

21.62 

1.6339 


5x5 

0.7970 

20.30 

-1.2266 

22.66 

2.2990 


6x6 

0.8002 

19.98 

-1.2260 

22.60 

3.0547 

6 

lxl 


>100 


>100 

0.9212 


2x2 

0.9357 

6.43 

-1.1711 

17.11 

1.1272 


3x3 

1.0349 

3.49 

-0.9648 

3.52 

1.4 82 8 


4x4 

1.0236 

2.86 

-0.9564 

4.36 

1.9553 


5x5 

1.0278 

2.78 

-0.9667 

3.33 

2.5755 


6x6 

1 . 0284 

2.84 

-0.9653 

3.47 

3.3316 

7 

lxl 


>100 


>100 

1.0463 


2x2 

0.9252 

7.48 

-1.2021 

20.21 

1.2503 


3x3 

1.0422 

4.22 

-0.9718 

2.82 

1.5904 


4x4 

1.0185 

1.85 

-0.9836 

1.64 

2.0833 


5x5 

1 .0176 

1.76 

-0.9794 

2.06 

2.6995 

12 

6x6 

1.0200 

2 . 00 

-0.9791 

2.09 

3.4563 


lxl 


>100 


>100 

1.9069 


2x2 

0 . 9057 

9.43 

-1.2149 

21.49 

2.1113 


3x3 

0.9799 

2.01 

-1.0127 

1.27 

2.4536 


4x4 

0.9926 

0.74 

-1.0061 

0.61 

2.9403 


5x5 

0 . 9962 

0.38 

-1.0021 

0.21 

3.5512 


6x6 

0.9986 

0 . 14 

-1.0017 

0.17 

4.3225 


NRI : Number of Gaussian quadrature points for regular integral 


NS I : Number of Gaussian quadrature points for singular integral 


Exact solution : u=l.E-06, t=-1.00 


Table 4.3 Beam in Uniform Tension, 6-node Triangular Element: Comparison 
of Error for Displacement and Traction, and CPU Time 


74 





















NRI 

NS I 

u A (xiCf 6 ) 

Error (%) 

C B 

Error (%) 

CPU Time 
( second ) 

lxl 

lxl 


>100 


>100 

0.1926 


2x2 


>100 

-0.3553 

64.47 

0.4299 


3x3 


>100 

-0.3430 

65.70 

0.8292 


4x4 


>100 

-0.3441 

65.59 

1.3935 


5x5 


>100 

-0 .3441 

65.59 

2.1057 


6x6 


>100 

-0 . 3441 

65.59 

2.9901 

2x2 

lxl 

1.1266 

12.66 

-0.6751 

32.49 

0.2595 


2x2 

0.8189 

18.11 

-1.2086 

20.86 

0.4985 


3x3 

0 . 8074 

19.26 

-1.2429 

24.29 

0.3954 


4x4 

0 . 8086 

19.14 

-1.2390 

23.90 

1.4562 


5x5 

0 . 8086 

19.14 

-1.1979 

19.79 

2.1756 


6x6 

0 . 8086 

19.14 

-1.2497 

24.97 

3.0527 

3x3 

lxl 

0 . 9108 

8 . 92 

-1 . 3704 

37.04 

C . 37 15 


2x2 

1.0389 

3.89 

-0 . 9424 

5.76 

0.6111 


3x3 

1.0228 

2.28 

-0 . 9710 

2.90 

1.0103 


4x4 

1.0210 

2 .10 

-0 . 9736 

2.64 

1.5556 


5x5 

0 . 9814 

1 .86 

-0 . 9738 

2.62 

2.2889 


6x6 

1.0208 

2.08 

-0.9738 

2.62 

3.1837 

4x4 

lxl 

0 . 9283 

7 .17 

-1.3788 

37.88 

0.5236 


2x2 

1 . 0187 

1 . 87 

-0 . 9648 

3.52 

0.7674 


3x3 

1 .0028 

0.28 

-0.9959 

0.41 

1.1598 


4x4 

0 . 9984 

0.16 

-0 . 9984 

0.16 

1.7192 


5x5 

0 . 9984 

0 .16 

-0 . 9990 

0.10 

2 .4486 


6x6 

0 . 9984 

0.16 

-0 . 9988 

0.12 

3.3228 

5x5 

lxl 

0 . 9278 

7.22 

-1.3453 

34.53 

0.7288 


2x2 

1 . 0139 

1.39 

-0 . 9666 

3.34 

0.9667 


3x3 

1 . 0014 

0 . 14 

-0.9971 

0.29 

1.3594 


4x4 

0 . 9997 

0 . 03 

-1.0003 

0.03 

1.9274 


5x5 

0 . 9997 

0 .03 

-1.0006 

0.06 

2.6425 


6x6 

0 . 9997 

0 .03 

-1.0006 

0.06 

3.4846 

6x6 

lxl 

0 . 9275 

7.25 

-1.3459 

34.59 

0.9826 


2x2 

1.0143 

1.43 

-0.9660 

3.40 

1.2132 


3x3 

1 . 0018 

0.18 

-0.9969 

0.31 

1.6159 


4x4 

1.0002 

0 . 02 

-0 . 9997 

0.03 

2.1875 


5x5 

1.0000 

0 . 00 

- 1.0000 

0.00 

2.8705 


6x6 

1.0000 

0 . 00 

- 1.0000 

0.00 

3.7663 

NRI : Number of Gaussian quadrature points for regular integral 


NS I : Number of Gaussian quadrature points for singular integral 
Exact solution : u=l.E-06, t=-1.00 


Table 4.4 Beam in Uniform Tension, 8-node Quadrilateral Element: Comparison 
of Error for Displacement and Traction, and CPU Time 
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NRI 

NSI 

u A ( * 10" 6 ) 

3 

3x3 

0 . 9426 


4x4 

0 . 9438 


5x5 

0 . 9438 


6x6 

0 . 9439 

4 

3x3 

0 . 9260 


4x4 

0 . 9272 


5x5 

0 . 9273 


6x6 

0 . 9273 

6 

3x3 

1.0034 


4x4 

1 . 0029 


5x5 

1 . 0030 


6x6 

1 . 0030 

7 

3x3 

0 . 9992 


4x4 

1 . 0028 


5x5 

1 .0029 


6x6 

1 .0028 

12 

3x3 

0 . 9974 


4x4 

1 . 0000 


5x5 

1.0000 


6x6 

1 . 0000 

13 

3x3 

0 . 9974 


4x4 

1 . 0000 


5x5 

1.0000 


6x6 

1 .0000 


Error (%) 


5.74 

-1.0365 

5.62 

-1.0418 

5.62 

-1.0415 

5 . 61 

-1.0418 

7 . 40 

-1.0529 

7.28 

-1.0579 

7.27 

-1.0582 

7 . 27 

-1.0582 

0.34 

-0.9928 

0.29 

-0.9978 

0.30 

-0.9981 

0 . 30 

-0.9981 

0.08 

-0.9923 

0.28 

-0.9972 

0.29 

-0.9975 

0.23 

-0.9975 

0.26 

-0.9948 

0.00 

-0.9997 

0.00 

-1.0000 

0.00 

-1.0000 

0.26 

-0.9948 

0.00 

-1.0000 

0.00 

-1.0000 

0.00 

-1.0000 


Error (%) 

CPU Time 
( second ) 

3 . 65 

0.6251 

4.18 

1.2028 

4.15 

0.8792 

4.18 

1.6006 

5.29 

0 . 6661 

5.79 

0 . 9229 

5.82 

1.2458 

5.82 

1 . 6534 

0.72 

0.7982 

0.22 

1.0496 

0.19 

1.3701 

0 .19 

1.7863 

0.77 

0.8463 

0.28 

1.0990 

0.25 

1.4222 

0.25 

1.8274 

0.52 

1.2506 

0.03 

1.5061 

0.00 

1.8256 

0.00 

2.2321 

0.52 

1.3041 

0.00 

1.5561 

0.00 

1.8756 

0.00 

2.2818 


NRI : Number of Gaussian quadrature points for regular integral 
NS I : Number of Gaussian quadrature points for singular integral 
Exact solution : u=l.E-06, t=-1.00 


Table 4.5 Beam in Uniform Tension, 3-node Triangular Element, Double Region: 
Comparison of Error for Displacement and Traction, and CPU Tune 
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u^( x 10' 6 ) 

Error (%) 


Error (%) 

CPU Time 
( second ) 

3x3 

3x3 

1 

0.17 

- 1.0000 

0.00 

0.8452 


4x4 

imERfl 

0 . 17 

- 1.0000 

0.00 

1.3190 


5x5 

1.0017 

0.17 

- 1.0000 

0.00 

1.9215 


6x6 

1.0017 

0 . 17 

- 1.0000 

0.00 

2.6509 

4x4 

3x3 

1.0000 

0 . 00 

- 1.0000 

0.00 

0.9502 


4x4 

1 .0000 

0.00 

- 1.0000 

0.00 

1.4150 


5x5 

1 .0000 

0 . 00 

- 1.0000 

0.00 

2.0153 


6x6 

1 . 0000 

0 . 00 

- 1.0000 

0.00 

2.7682 

5x5 

3x3 

1 . 0000 

0.00 

- 1.0000 

0.00 

1.0912 


4x4 

1.0000 

0.00 

- 1.0000 

0.00 

1.5542 


5x5 

1.0000 

0 . 00 

- 1.0000 

0.00 

2.1523 


6x6 

1 .0000 

0 .00 

- 1.0000 

0.00 

2.8998 

6x6 

3x3 

1.0000 

0.00 

- 1.0000 

0.00 

1.2692 


4x4 

1.0000 

0.00 

- 1.0000 

0.00 

1.7257 


5x5 

1.0000 

0.00 

- 1.0000 

0.00 

2.3466 


6x6 

1 . 0000 

0.00 

- 1.0000 

0.00 

3.0907 


NRI : Number of Gaussian quadrature points for regular integral 
NS I : Number of Gaussian quadrature points for singular integral 
Exact solution : u=l.E-06, t=-1.00 


Table 4.6 Beam in Uniform Tension, 4-node Quadrilateral Element, Double Region: 
Comparison of Error for Displacement and Traction, and CPU Time 
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NRI 

NS I 

u A ( X 10' 6 ) 

Error (%) 


Error (%) 

CPU Time 
( second ) 

4 

3x3 

0.7582 

24.18 

-1.2355 

23.55 

2.3109 


4x4 

0.7391 

26.09 

-1.2142 

21.42 

3.2921 


5x5 

0 . 7420 

25 . 80 

-1.2183 

21.83 

4.5192 


6x6 

0 .7448 

25.52 

-1.2182 

21.82 

6.0571 

5 

3x3 

1 . 0412 

4 . 12 

-0 . 9650 

3.50 

2 . 8737 


4x4 

1 . 0299 

2 . 99 

-0.9536 

4.64 

3.8452 


5x5 

1 . 0304 

3.04 

-0.9646 

3.54 

5.0775 


6x6 

1 .0311 

3 .11 

-0.9638 

3.62 

5.6007 

7 

3x3 

1 .0458 

4.58 

-1.0546 

5.46 

3.1014 


4x4 

1 . 0212 

2.12 

-1.0357 

3.57 

4.0894 


5x5 

1.0236 

2.36 

-1 .0441 

4.41 

5.2954 


6x6 

1 . 0257 

2.57 

-1 .0443 

4.43 

6.8005 

12 

3x3 

1 .0197 

1 . 97 

-1.0087 

0.87 

4.8588 


4x4 

0 . 9932 

0 . 68 

-0.9881 

1.19 

5.8394 


5x5 

0 . 9968 

0 . 32 

-1.0017 

0.17 

7.0863 


6x6 

0 . 9991 

0 .09 

-0.9994 

0.06 

3.5758 

13 

3x3 

1 .0192 

1 . 92 

-1.0086 

0.86 

5.1563 


4x4 

0 . 9928 

0.72 

-0.9883 

1.17 

5.1327 


5x5 

0 . 9964 

0.36 

-0.9948 

0.52 

7.3690 


6x6 

0 . 9987 

0 . 13 

-0.9984 

0.16 

3.8721 


NRI : Number of Gaussian quadrature points for regular integral 
NS I : Number of Gaussian quadrature points for singular integral 
Exact solution : u=l.E-06, t=- 1 .00 


Table 4.7 Beam in Lniform Tension, 6-node Triangular Element, Double Region: 
Comparison of Error for Displacement and Traction, and CPU Time 
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■6 \ 



Error (%) 

CPU Time 

NR I 

NS I 

u ( * 10 ) 

A 

Error (%) 



3x3 

3x3 

1 .0179 

1.79 

-0 . 9653 

3.47 

1.9954 


4x4 

1.0170 

1.70 

-0 . 9680 

3.20 

3 . 1194 


5x5 

1 . 0170 

1.70 

-0 . 9683 

3.17 

4.5595 


6x6 

1 .0169 

1.69 

-0 . 9683 

3.17 

6.3246 

4x4 

3x3 

1.0017 

1.70 

-0 . 9955 

0.45 

2.3013 

4x4 

1 . 0012 

0 . 12 

-0 . 9986 

0.14 

3.4235 


5x5 

1.0012 

0 . 12 

-0 . 9988 

0.12 

4.8600 


6x6 

1 . 0012 

0.12 

-0 . 9989 

0.11 

6.6606 

5x5 

3x3 

1.0006 

0.12 

-0.9973 

0.27 

2.7091 

4x4 

0 . 9999 

0.01 

-1 . 0004 

0.04 

3.8227 


5x5 

0 . 9996 

0.04 

-1.0007 

0.07 

5.2572 


6x6 

0 . 9996 

0.04 

-1.0007 

0.07 

7.0385 

6x6 

3x3 

0 . 9993 

0 . 07 

-0 . 9966 

0.34 

3.2110 


4x4 

1.0000 

0.00 

-0.9997 

0.03 

4 . 3337 


5x5 

1 . 0000 

0.00 

-1 . 0000 

0.00 

5.7791 


6x6 

1 . 0000 

0 . 00 

- 1.0000 

0.00 

7 . 5443 

NRI : Number of Gaussian quadrature points for regular integral 


NS I : Number of Gaussian quadrature points for singular integral 

Exact solution : u=l.E-06, t=-1.00 





Table 4.8 Beam in Uniform Tension, 8-node Quadrilateral Element, Double Region: 
Comparison of Error for Displacement and Traction, and CPU Tune 
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4.2 Beam Subjected to Pure Bending 


The problem of a beam in pure bending was solved to investigate the 
convergence of each element as a function of number of elements used. The 
geometry and boundary conditions are shown in Figure 4.6. For each element, 
three different meshes were used as shown in Figures 4.7 to 4.10. The exact 
solution for the displacements are [31] 


_ vM 


u * = ~eT xy 


Uy = ' EI ( z2 ' Vx2+ Vy2 ) 


_U 
EI 


u z=T?7 yz 


(4.2) 


and the traction at point B in Figure 4.6 is 


h = a (4.3) 

where M is the applied moment and I is the moment of inertia with respect to the z 
axis. The displacements at point A and the traction in direction z at point B are 
compared with these exact solutions. The results are listed in Tables 4.9 to 4.12. 
The 3-node triangular element and 4-node quadrilateral element do not converge to 
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Figure 4.6 Geometry and Boundary Conditions for Beam in Pure Bending 
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Figure 4.7 3-Node Triangular Element Meshes for Beam in Pure Bending 
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Figure 4.8 4-Node Quadrilateral Element Meshes for Beam in Pure Bending 
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Exact solution : u x = 1 .OOOOOE-07. u = 1.06667E-05, u , = -2.66667E-06, t =1.00 

y z z 

Table 4.10 Beam in Pure Bending, 4-node Quadrilateral Element.Comparison of Error for Displacementsand Traction 






Figure 4.9 6-Node Triangular Element Meshes for Beam in Pure Bending 
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Table 4.1 1 Beam in Pure Bending, 6-node Triangular Element Comparison of Error for Displacements and Traction, and 
CPU Time 
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Table 4.12 Beam in Pure Bending, 8-node Quadrilateral Element Comparison of Error for Displacements and Traction, 
and CPU Time 




the exact solution even for the highest order of Gaussian quadrature and for the 
finest mesh. These results lead to the conclusion that the linear elements should 
not be used for problems in which the order of traction and displacement 
distributions are higher than linear. The 6-node triangular element and 8-node 
quadrilateral element both converge to the exact solutions either by increasing the 
number of Gaussian quadrature points or by refining the mesh. It is also concluded 
that the 8-node element is superior to the 6-node element In order to obtain results 
with less than 1% error in both the traction and displacement, the 6-node element 
requires at least 17.9984 seconds using the 7 point Gaussian quadrature for the 
regular integral and 6x6 for the singular integral on the 42 node mesh but the 8- 
node element only needs 7.5160 seconds by choosing 4x4 Gaussian quadrature for 
the regular integral and 3x3 for the singular integral on the 56 node mesh. 

The multi-domain technique was also applied using the 8-node element. 
Two different meshes were tested, as shown in Figure 4.11. The results in Table 
4.13 show that the double region mesh can achieve an equivalent accuracy by 
consumming less CPU time than that the single domain mesh does. This is because 
the system matrices for the two subregions are, in fact, separated and thus for the 
number of nodes the double region mesh needs less calculation than the single 
domain mesh. 
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Interface 



(a) Total no. of nodes: 40 ( 20 nodes for each subregion) 
Total no. of Elements: 12 (6 elements for each subregion) 


Interface 



(b) Total no. of nodes: 64 ( 32 nodes for each subregion) 
Total no. of Elements: 20 (10 elements for each subregion) 


Figure 4. 1 1 8-Node Quadrilateral Element, Two-Subregion Meshes for Beam 
in Pure Bending 
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4.3 Circular Buried Crack Under Uniform Tension 


The stress intensity factor of a circular crack buried in an infinite body 
subjected to far field uniform distributed traction is studied in this section. The 
exact solution is [32] 


K r -W# 


(4.4) 


where a is the applied stress and a is the crack radius. For this problem Kjj and 
K ttt are zero since the load is perpendicular to the crack surface. Six different 
meshes with either different number of nodes or different boundary conditions were 
studied. The overall dimensions of these meshes were chosen large enough to 
simulate the infinite medium. 

The first mesh is a two-subregion model which describes the whole domain 
of the problem. Each subregion consists of 52 elements and 150 nodes in which 88 
nodes belong to the interface which bonds the two subregions as shown in Figure 
4.12. The crack face which is kept traction free. Each crack face is formed by 
eight 8-node quadrilateral quarter point elements and eight 6-node triangular 
elements. Traction singular elements are placed along the crack front The uniform 
distributed tractions o in direction z are applied on the plane parallel to the interface 
plane of one subregion mesh, and the same plane for another subregion is then 
fixed in the z-direction with node A fixed in direction x and node B fixed in 
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C -L 




Figure 4.12 Element Mesh of Double Region Model for Buried Circular Crack 
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Error (%) 



Fi§-4.13 Average Error in Computed KI by Different L/A for Varied Poisson Ratio 
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direction y, as showm in Fig.4.12, to prevent rigid body translation and rotation. 
Ingraffea and Mann [33] pointed out that the ratio of the length L of the crack front 
element to the crack length a affects the accuracy of the results. Therefore, 
convergence for different L/a ratios under varied Poisson's ratio v was studied first. 
Five different ratios of L/a varying from 0.1 to 0.5 were applied for Poisson's 
ratios varying from 0.0 to 0.4. The results are shown in Figure 4.13 revealing that 
the absolute error is confined to within 5% when the ratio of L/a is in between 0.26 
and 0.34 for Poisson's ratios ranging from 0.1 to 0.4. Thus, the ratio of L/a of the 
models used in this paper is chosen to be 0.3. 

In the second example the problem is modeled using the coarse double 
region mesh shown in Figure 4.14. It consists of 73 nodes and 25 elements. The 
ratio of L/a is 0.3. The maximum error in the calculated Kj was found to be 
0 . 23 %. The CPU time for this case was 31 seconds. 

The third and fourth meshes used are the same as the ones above except that 
only half of the subregion is taken into consideration. The nodes belonging to the 
interface are now constrained in the z-direction as shown in Figure 4. 15 and 4.16 to 
simulate the correct symmetric boundary conditions. The calculated Kj is almost 
constant along the crack front and the absolute errors were 1.6% and 0.27% for the 
mesh of Figure 4.15 and the mesh in Figure 4.16, respectively. The CPU times 
were 71 seconds and 20 seconds, respectively. 

The fifth mesh is the half cut of the third mesh with an additional cut plane 
y=0 which is fixed in the y-direction as shown in Figure 4.17. The 6-node 
triangular elements around the crack tip on this cut plane are converted to be traction 
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Figure 4.14 Element Mesh of Double Region Model for Buried Circular Crack 
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Figure 4. 15 (a) Element Mesh for Half Domain of Buried Circular Crack 
(b) Illustration of Elements around Crack 
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(b) 


No. of Elements: 25 



Crack Front 


Traction Singular 
Element 


Quarter Point 
Element 


Figure 4. 1 6 (a) Element Mesh for Half Domain of Buried Circular Crack 
(b) Illustration of Elements around Crack 
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Traction Singular 



Element 
Crack Front 

Quarter Point 
Element 


Total number of nodes: 202 
Total number of elements: 72 


Figure 4.17 Element Mesh of Quarter Domain for Buried Circular Crack 


99 




Figure 4.18 Element Mesh of One- Eighth Domain for Buried Circular Crack 
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Mesh 

(fig.) 

Mesh 

Type 


4 . 12 


Double 

Region 


4 .14 

Double 

Region 


4.15 


Single 

Region 


4 . 16 


Single 

Region 


4.17 

Single 

Region 


4.18 


Single 

Region 


No. of 
Node 


300 


146 


150 


73 


202 


110 


No. of 
Element 


104 


50 


52 


25 


72 


40 


Ki 


1 . 14604 


1.13101 


1.14647 


1.13141 


1.15534 


1.15659 


Error 

(%) 


1 . 6 


0.23 


1 . 6 


0.27 


2.4 


2.5 


CPU Time 
(Second) 


131 


31 


71 


20 


133 


41 


Exact solution : K t = 1.128379 

where crack radius a=l, applied stress <7=1. 


Table 4. 14 Circular Crack Buried in Infinite Body Under Uniform Tension: 
Comparison of Error for Kj and CPU Time for Different Meshes 
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singular elements by relocating their middle nodes on the two sides to the quarter 
position to capture the inverse square root singularity in stresses. This mesh is 
composed of 202 nodes and 72 elements. The differences between the analytic 
solution of Kj and the calculated results ranged from 1.2% to 2.3%. This mesh 
uses 133 seconds of CPU time. 


The final mesh models one-eighth of the domain. Three planes, x=0, y=0, 
and z=0, are fixed in the x, y, and z direction, respectively, to represent the correct 
boundary conditions. The maximum error in Kj is 2.5%. 

The results of this section are summarized in Table 4.14. From the table we 
also observe that the last two cases which model only part of the whole domain 
have the largest error. This is probably due to the asymmetry of the meshes. 
Suprisingly, the coarsest mesh (73 nodes) leads to the most accurate result. Since 
the Kj is constant along the crack front and the quadratic element on both side of the 
crack front can model the circular shape exactly, an accurate result can still be 
obtained by using only a few elements to model the crack. This also counts for the 
reason that the coarsest mesh can have such excellent results. 


4.4 Circular Buried Crack Inclined at 30 Degrees Under Uniform Tension 

This problem is used to ensure the ability of the Boundary Element Method 
to calculate stress intensity factors for mixed mode fracture problems. A circular 
crack deforms in three modes when the normal of the crack surface is not parallel to 
the direction of the applied load. The exact solutions for the stress intensity factors 
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Figure 4.19 Circular Crack Buried in Infinite Body Inclined at an Angle y With 
Respect to the Direction of Applied Load 
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are [33] 


K i= l|-asin 2 V7mJ 


K n 


k(2-v) 


(asiny cosy)/ia 


COSO) 


K m = 


^2-v) ( asin T cosy Via 


srnco 


(4.5) 


where a is the applied stress, a the crack radius, y the angle between the direction of 
applied load and the normal of the crack plane, and co is the angle as is shown in 
Figure 4.19. Only a double region mesh can model this problem appropriately 
since symmetry no longer exists. The mesh is composed of two parallelepiped 
shaped subregions as shown in Figure 4.20. Each subregion is made up of 52 
quadratic elements and 150 nodes including 88 nodes belonging to the interface. 
Quarter point elements and traction singular elements are used along the crack front 
as described in the previous section. The results are shown in Figures 4.21 to 4.23 
in which the stress intensity factors are normalized. The calculated Kj are almost 
constant along the crack front. The error ranges from 1.17% to 1.19%. The 
maximum error for Kjj and Kjjj are 2.34% and 5.71%, respectively. The total 
CPU time is 131 seconds. 
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Figure 4.20 Element Mesh of Double Region Model for Buried Circular Crack 
Inclined at 30 Degrees 
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Figure 4.22 Comparison of Calculated KII with Exact Solution for 
Circular Crack Inclined at 30 Degrees 
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4.5 Semi-circular Surface Crack Under Uniform Traction 


The purpose of this section is to check the accuracy of the method for 
surface crack problems. Since an exact solution for this problem does not exist, the 
results are compared with the work of Tada [33], and Newman and Raju [34], 
Tada presents the stress intensity factor for a semi-circular surface crack in a 
semi-infinite body as shown in Figure 4.24 as 


Kj(0) =|-aVia F(6) 


(4.6) 


where 


F(9) = 1.211 - 0.186 VsinQ (10'<9<170‘) 


and c is the applied stress, a the crack radius, and 0 is the angle measured from the 
surface as shown in Figure 4.24. Note that the Ki is not constant along the crack 
front. For the same problemNewman and Raju predict 


K x = cfi ta (rj=) F(9) 


(4.7) 
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where 


F(0) = 1.04 [ 1+0.1 (l-sin0) 2 ] 


and Q = 2.464 for a semi-circular surface crack. A double region model, as shown 
in Figure 4.25, is used in the problem. Each subregion consists of 72 quadratic 
elements and 202 nodes. 

The comparison of the results with the work of Tada and Raju and Newman 
is shown in Figure 4.26. The present results are quite consistent with these two 
predictions when the 9 is between 35’ and 90*. The max im u m difference appears 
at the surface with the error of about 5.7%. This is because the Kj is calculated 
assuming of plane strain conditions and an inverse square root singularity of 
stresses. At the points where the crack intersects the free surface these conditions 
do not apply. 


Ill 
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Figure 4.25 Element Mesh for Double Region Model of Surface Crack 
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Figure 4.26 Comparison of Stress Intensity Factors for a Circular 
Surface Crack in a Semi-infinite Body Under Uniform 
Tension 
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CHAPTER FIVE 


STRESS INTENSITY FACTOR ANALYSIS OF THE INNER 
RACEWAY OF HIGH SPEED BEARINGS 


In this chapter, the geometry of the engine bearing is described and the 
loadings on the inner raceway (including the hoop stresses and the Hertzian contact 
load) are calculated. The stress intensity factors for several semi-circular surface 
cracks of different lengths and inclinations are presented as functions of the location 
of the indenters. 

5. 1 Geometry and Applied Loading 

The bearing analyzed in this report is a high performance engine bearing 
which is used on the main shaft of an aircraft. The dimensions and the geometry 
of the bearing are shown in Table 5.1 and Fig.5.1, respectively. 

The bearing consists of 28 ball rollers. To simulate the passage of each ball 
only l/28th of the inner raceway is modeled. Since the radius of the inner raceway 
is large compared to the other dimensions of the part, the curvature is neglected and 
the mesh is modeled as a block with flat surface as shown in Fig. 5. 2. 

The external loadings considered in the analysis are the hoop stresses and 
the Hertzian load. The hoop stresses due to the rotation and the shrink fit are taken 
from [22] and are given by 
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Shaft 


Inner raceway 


Outer raceway 


Bearing length L 


Ball bearing 


No. of ball bearings 


Shaft speed 


M50 steel 
material properities 


Inner radius a =2.0 inch 
Outer radius b s = 2.30233 inch 
Inner radius a j =2.3 inch 
Outer radius b j = 2.6 inch 
Inner radius a 0 =3.1 inch 
Outer radius b 0 = 3.35 inch 
0.57322 inch 


Radius R = 0.25 inch 


28 


25,500 rpm 
E = 3 .OX 10 7 psi 
p = 0.288 lb/in 3 


v = 0.3 

IU= 18 ksWin 


Table 5.1 Dimensions and Material Properities of Typical Ball Bearing for 
Aircraft Engines 
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Figure 5. 1 Geometry of Bearing 
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(5.1) 


a 


ee 



l+3v 

3+v 


r 2 ] + P 


(b./r ) 2 + 1 
(b/a.) 2 - 1 


where co is the angular velocity of the shaft and P is the pressure existing between 
the shaft and the inner raceway which is equal to 


E5 n (a? - a 2 ) ( b. - a 2 ) 
a i 2a 2 (bf - a? ) 


(5.2) 


where 8 n is the difference between the outer radius of the shaft and the inner radius 
of the inner raceway during rotation at speed co. It can be obtained by 


S n =S-AS (5-3) 

where 5 is the original shrink fit at 0 rpm which is equal to 0.00233 inch in this 
case and AS is the difference in the radial displacement between the inner radius of 
the inner raceway, u^, and the outer radius of the shaft, u r s , under the rotation at 
speed co, i.e.. 
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A5 = uj - IV s 

= P© 2 ^ 
4E 


[(3+v)bj +(l-v)af]- 


pco 2 b s 

4E 


[ (3+v)b s 2 + (1 - v)aj] (5.4) 


Applying the data of Table 5.1 to the above equations results in 


<j 09 = 1973.61[ 3.5986 + ] 


r 2 

CD 2 x 10' 5 [ 3.6290 - 0. 1941(r 2 ) + 2J2*! j ( ps i) 

r 2 


(5.5) 


where 0) is in rpm. The maximum hoop stress of the inner raceway is at the inner 
radius. The hoop stresses range from 38 ksi to 44 ksi when the rotation speed is 
25,500 rpm. 

The contact region between a sphere and a flat surface is circular and the 
stress distribution is [30] 



(5.6) 


where a is the radius of the contact area and p 0 is the maximum stress which can be 
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calculated from 


Po 


3 P T 
2 7ta 2 


(5.7) 


where Pj is the total load. The radius is 


a = 


3 /3P t R(1-v^) 
V 2E 


(5.8) 


where R is the radius of the sphere. As discussed in [6] the experimental data 
indicates that for the bearing considered in the present analysis the maximum 
Hertzian stress is 285 ksi. Subsequent calculations will be based on this value. 

The distribution of the Hertzian contact loading is distorted when the 
spherical body is rolling over a surface containing a surface breaking crack [31]. 
However, the influence of the subsurface crack on the Hertzian stress distribution 
has been shown [31] to be insignificant. The interaction effects of the subsurface 
crack on the Hertzian distribution are thus ignored. 

Since the Hertzian stress fields can be solved in closed form, superposition 
is applied in order to decrease the number of nodes needed to model the Hertzian 
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contact load accurately. The superposition method is illustrated in Fig.5.3. The 
stress intensity factors of a surface crack subjected to the Hertzian load are equal to 
the stress intensity factors of the same crack loaded with the negative of the stresses 
produced by the Hertzian loading. The stresses produced in the interior of the plate 
by the hertzian contact loading are calculated by integrating the stresses due to a 
concentrated force acting on the boundary of a semi-infinite solid [22] as shown in 
Fig.5.4. The hoop stresses are directly applied on the the end of the plate. The 
total stress intensity factor is then the superposition of the stress intensity factor due 
to the hoop stresses and the one resulting from the Hertzian load. It should be 
noted that the depth of the plate is assumed to be large enough so that the 
distribution of loading remains Hertzian. 


5.2 Stress Intensity Factor of Circular Subsurface Crack in the Inner Raceway of 
the Engine Bearing 

This section presents the results of a quasi-static stress intensity factor 
analysis of a typical ball bearing which is used as a support for the main shaft of 
aircraft engines. All the calculations are based on the model shown in Fig.5.2 
which consists of an axial hoop stress and a Hertzian load ( contact radius a and 
maximum intensity p 0 ) interacting with a semi-circular crack (radius, or length 1 ) 
inclined at an angle <(>. The dynamic effect is ignored and the distance x, between 
the center of the Hertzian distribution and the crack mouth is changed incrementally 
to simulate the passage of each ball bearing. 

The stress intensity factors for this problem vary with position along the 
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Figure 5.3 Illustration of Superposition Method 
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Figure 5.4 Concentrated Force acting on the Boundary of a Semi-infinite Body 
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crack front (9 ). A typical variation of the mode I stress intensity factor at 9 = 90° 
with roller position and crack length for a vertical crack ( <}> = 0° ) is shown in 
Fig.5.5 for p Q = 285 ksi ( a = 0.00679 inch ). When the roller is at a distance 
greater than four times the contact length from the crack mouth, the stress intensity 
factor is a constant which results from the axial stress. As the roller gets closer, the 
compressive stress arising from the Hertzian load decrease the Kj stress intensity 
factor. When the load is on the crack mouth the stress intensity factors becomes 
negative for the cracks which are shorter than 0.005 inch. The negative value of Kj 
indicates the closure of the crack. For longer cracks, the decrease of Kj diminishes 
since the crack tip is beyond the range of the highly compressive Hertzian stress 
field. 


The variation of the mode II stress intensity factor for the same loading 
condition is shown in Fig.5.6. The value of Kn is zero when the roller is far from 
the crack mouth. As the roller approaches the crack, Kjj starts to increase and 
reaches a maximum value when the load reaches the crack mouth. As the roller 
crosses to the other side of the crack, Kjj abruptly changes sign and decreases to a 
minimum value equal in magnitude to the previous maximum. For small cracks this 
change is very abrupt, but for large cracks the change is more gradually. As 
pointed out in [6], these abrupt variations in Kq and Kjj may significantly affect the 
propagation of short cracks. 

The stress intensity factors for 9 = 45° are shown in Fig.5.7 to 5.9. It is 
observed that the variations of Kj and Kn did not differ significantly from 9 = 90°. 
Another stress intensity factor Km is observed based on the local coordinate 
system moving along the crack front. The Km behaves as Kn- Fig.5.10 to 5.12 
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KI (ksiVin) 


Roller Position (x/a) 

Figure 5.5 Variation of KI at 6 = 90° as a Function of Roller 
Position with Different Crack Length 
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Roller position (x/a) 

Figure 5.6 Variation of KII at 6 = 90° as a Function of Roller 
Position with Different Crack Length 
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KI (ksWin ) 


Roller position 

Figure 5.7 Variation of KH at 9 45’ as a Function of Roller Position with 
Different Crack Length 
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KII (ksiVin ) 



Roller position (x/a) 

Figure 5.8 Variation of KH at 0 =45° as a Function of Roller Position with 
Different Crack Length 
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Kill (ksiVin ) 



Figure 5.9 Variation of Kill at 0 = 45° as a Function of Roller 
Position with Different Crack Length 
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show the stress intensity factors for different locations along the crack front. It can 
be seen from Fig.5. 10 that the magnitude of Kj does not change significantly. The 
magnitudes of Kn and Km change along the crack front but however it can be 
calculated from Fig.5. 1 1 and 5. 12 that the square root of the sum of the squares of 
Kq and Km are almost constant along the crack front. A comparison of the results 
for <(> = 0° obtained in the presented analysis with those presented in Mendelson and 
Ghosn [6] revealed that the magnitudes and variations in Kj are similar. However, 
the mode II stress intensity factors differ significantly. As seen in Fig.5. 5, for 1 = 
0.02 inch the maximum value of Kjj is approximately 1.5 ksi Vin, while in 
reference [6] it is approximately 10 ksiVin. 

Fig.5. 13 and 5.14 present the stress intensity factors of K{ and Kji at 
9=90° , respectively, for cracks inclined at <|> = 30° for p 0 = 285 ksi. An increase of 
Ki is observed when the roller passes to the right hand side of the crack mouth for 
short cracks. This is because the Hertzian load causes the inclined crack surfaces 
apart when it is passing over the crack mouth. The value of Kq before the roller 
crossing over the crack from the left is much greater than the Kjj after the roller 
moving to the right of the crack mouth since the Hertzian load is pushing the left 
crack surface sliding along the right crack surface when the roller is on the right 
hand side of the crack. 


Fig.5. 15 and 5.16 show Kj and Kjj variations for several inclinations of a 
crack length 1 = 0.02 inch. High values of Kjj is observed for <f» = 15° and <(> = 30°. 
The mode I fracture toughness of M50 steel, which is used for this type of bearing, 
is of the order of 18 ksWin. Assuming the mode II toughness is of comparable 
magnitude, these results indicate that this applied loading would lead to large 
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propagation rates for cracks inclined at 30°. 
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Roller position (x/a) 

Figure 5.10 Variation of KI as a Function of Roller Position with Different 
Angle 0 
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Figure 5.11 Variation of Kn as a Function of Roller Position with Different 
Angle 9 
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Angle ® 
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Roller Position (x/a) 

Figure 5.13 Variation of KI at 9 = 90° as a Function of Roller 
Position with Different Crack Length 
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14 Variation of KII at 6 =’90° as a Function of Roller 
Position with Different Crack Length 
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Kl (ksiVin ) 



Roller position (x/a) 

Figure 5.15 Variation of KI at 9 = 90 ° as a Function of Roller 
Position with Different Crack Length 
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Roller position (x/a) 

Figure 5.16 Variation of KII at 9 = 90° as a Function of Roller 
Position with Different Inclined Angle 0 
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CHAPTER SIX 
CONCLUSIONS 


A preliminary stress intensity factor analysis of a typical high speed bearing 
was conducted using the Boundary Element Method. The results obtained in the 
present three-dimensional analysis suggest lower mode Kn stress intensity factors 
than those predicted by the two-dimensional analysis in [6]. This is due to the fact 
that the total load needed to produce the experimentally measured 285 ksi Hertzian 
stress using a semi-spherical contact area is much lower than that using a cylindrical 
contact (27.5 lbs instead of 1500 lbs ). This may be the reason why the predictions 
of Mendelson and Ghosns’ [6] analyses are overly conservative. High Kj and Kn 
values were observed for cracks inclined at 30°. These results indicate that the 
interaction of the Hertzian load would lead to large propagation rates for cracks 
inclined at 30°. 

Although the stress intensity factor data obtained from the analysis has not 
been reduced to a form suitable for life prediction, these preliminary results can 
provide a better understanding of the complex interactions between a surface crack, 
a moving Hertzian load, and an axial stress. 

As for the further work, an incremental crack growth analysis of elliptical 
cracks using the Boundary Element Method and a fatigue crack growth law would 
be the next step. Also, more factors which affect the stress intensity factors of a 
crack such as friction between the roller and the raceway, the dynamic effect, etc. 
could be taken into consideration. 
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APPENDIX A 

SHAPE FUNCTIONS FOR ISOPARAMETRIC ELEMENTS 


The shape functions for different elements are listed below corresponding to 
the elements shown in Fig.A.l. 

(1) 3-Node Linear Triangular Element: 


N1-S1 
N 2 = *2 
N 3 = ^3 

where * 3 = 1 

(2) 6-Node Quadratic Triangular Element: 


Nl-5l< 2*1-0 
N 2 = 5 2 (2*2-D 

N 3 = 4 3 (2^3-1) 
N 4 = 44^ 

N 5 = 4^3 
N 6 = 4*3*l 

where ?j + *3 = 1 
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(3) 4— Node Linear Quadrilateral Element: 


N 1 =(l-^ 1 )(l-^ 2 )/4 

n 2 = ( 1 +5i)(l ~%2 ) /4 
N 3 = (1+§i)( 1+§ 2 )/4 
N 4 = d- ^)(l + ^ 2 )/4 


(4) 8-Node Quadratic Quadrilateral Element: 

N 1 = '(l+^l+^X 1 ' ^l)( 1 - ^2 ) /4 
N 2 = -(l+^-^ 2 )( 1 + ^i)(l - £ 2 ) / 4 
N 3 = (-l+^ 2 )( 1 + ^)( 1 + % 2 ) / 4 

N 4 =-(1+5H 2 X 1 - ^i)(1+^ 2 )/4 

N5 = (l-^2xi-4 2)/ 2 

n 6 = ( 1 + ^iX1-42 2 V2 
N 7 = (l-^ 2 )(l+5 2 )/2 
N 8 = (l- 5 i)U-5 2 2 )/2 
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2 


*3 = '-« r «2 


2 



(1) 3-Node Linear 
Triangular Element 



(2) 6- Node Quadratic 
Triangular Element 



(3) 4- Node Linear 

Quadrilateral Element 





(4) 8-Node Quadratic 
Quadrilateral Element 


Figure A. 1 Isoparametric Elements 
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APPENDIX B 


TRANSFORMATION FUNCTIONS 


The integral on a boundary surface element with domain T can be expressed 


as 



r 


In order to overcome the 1/r singularity of the kernel the integral is first transformed 
from the Cartesian coordinate system to the parametric ^-coordinate system by 
using the shape functions 


“m 

* = 4 2 > x i 

i=l 

n m 

y = Z N i«i^) Y i 

i=l 


corresponding with the Jacobian 
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dx dx 
*1 ’ *2 


dy dy 


such that the integral becomes 



The integral on the isoparametric element is then divided into several triangles 
according to the location of the singular node. The singular integral on each triangle 
is carried out by using a polar coordinate system with its origin at the singular node 
such that 


^ = gi(r,6) 
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J 2 (r,0) = 




dr 

’ ae 

*2 

*2 

dr 

’ ae 


and 


"t 

-SJ K(r,9)J 1 (r,e)J 2 (r,e)drde 

i=l j 

i 


where nj is the total number of triangles in which the isoparametric element is 
divided and Tj is the corresponding domain. The Jacobian — r can remove the 
1/r singularity. In order to accomplish the integral numerically by the Gaussian 
quadrature the polar coordinate is again transfered to a system with both coordinates 
ranging from -1 to 1 by the transforming functions 


r = h.(r , 0 ) 
e = l.(r,0) 


148 



dr 3r_ 

dr ’ 30 


J 


3 


30 30 

3r ’ 30 


The integral then becomes 


"t 

-2J K(r,6)J 1 (r,0)J 2 (r,0)J 3 (r, 0 )drd0 

i=l j 


-2 


22 w^w^kO*,© 1 ^^.© )J 2 (^,0 )J 3 (^,0 )drd0 


i=l 


L a=l b=l 


where n a and n^ are the order of the Gaussian quadrature. The transforming 
functions for different element are illustrated in the following pages. 
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Figure B.l.a Illustration of Transformation of Coordinate System 
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Figure B.l.b Illustration of Transformation of Coordinate System 
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Figure B.l.c Illustration of Transformation of Coordinate System 
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Figure B.l.d Illustration of Transformation of Coordinate System 
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Figure B.l.e Illustration of Transformation of Coordinate System 
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Figure B. l.f Illustration of Transformation of Coordinate System 
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Figure B.l.g Illustration of Transformation of Coordinate System 
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Figure B.l.h Illustration 


E,l= rcos 0 +1/2 
^2 = r sin 9 +1/2 

J 2 = r 

e = — J-( e + 6) 

(r+l) 

1 4 cos (3ti/2-0 ) 

J = 2 

3 16cos(37t/2-0 ) 

-l< r_< 1 

-i < a.< i 


of Coordinate System 
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